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1.  Summary 


The  United  States  Air  Force  (USAF)  is  facing  a  potentially  large  wave  of  turbine  engine 
disc  replacement  costs  over  the  next  eight  to  ten  years  that  are  inconsistent  with  anticipated 
budgets.  Consequently,  the  Engine  Rotor  Life  Extension  (ERLE)  program  was  conceived 
by  the  AFRL  as  a  sound  science  and  technology  investment  that  offers  the  potential  for 
significant  cost-avoidance  by  extending  the  life  of  certain  life-limiting  components.  The 
concept  is  to  extend  the  life  of  these  components  by  recovering  the  conservatism  believed 
to  exist  in  current  design  and  life  management  practices,  without  increasing  risk,  by 
systematically  improving  and  more  effectively  integrating  a  number  of  life  management 
technologies — life  prediction,  non-destructive  inspection,  engine  health  monitoring, 
maintenance  and  repair.  Enhancements  in  engine  life  management  technology  would  also 
be  applicable  to  developmental  and  future  military  engines,  commercial  engines,  and  land- 
based  combustion  turbines  where  safety,  reliability,  and  cost  of  ownership  are  of 
paramount  importance.  SwRI  had  lead  a  team  with  unique  capabilities  to  enhance,  as  well 
as  integrate,  several  of  the  above  life  management  technologies.  Other  team  members 
include  Smiths  Aerospace,  The  University  of  Texas  at  San  Antonio,  and  Mustard  Seed 
Software. 

The  general  approach  and  technical  objectives  of  this  program  were  to:  1)  develop  and 
validate  a  new  family  of  physically-based,  deterministic  life  prediction  models  for  treating 
total  fatigue  life  including  crack  nucleation,  microcrack  growth  and  large  crack  growth; 

2)  develop  and  demonstrate  an  efficient  probabilistic  life  prediction  methodology  based  on 
the  stochastic  nature  of  each  of  the  above  phases  of  fatigue  life;  and  3)  explore  methods  for 
enhanced  engine  life  management  based  on  hybridisation  of  state-of-the-art  probabilistic 
life  prediction  and  classical  engine  health  monitoring. 

The  program  has  developed  and  validated  a  new  deterministic  fatigue  crack  nucleation  and 
growth  (FaNG)  model.  Model  predictions  were  found  to  be  in  good  agreement  with 
measured  fatigue  lives  for  notched  Ti-6A1-4V  specimens  over  a  wide  range  of  loading 
conditions.  The  effectiveness  of  the  FaNG  model  has  been  demonstrated  by  predicting  the 
beneficial  effect  of  residual  stresses  introduced  by  shot  peening  and  low  plasticity 
burnishing  on  enhancing  fatigue  life.  Methods  have  also  been  developed  and  implemented 
to  perform  probabilistic  sensitivity  analyses  to  rank  the  impact  of  key  random  variables  on 
component  reliability.  Probabilistic  tools  were  also  developed  to  automatically  identify 
missions  from  engine  usage  data  from  aircraft  flight  data  recorders,  as  well  as  use  this 
information  in  combination  with  the  newly  developed  fatigue  models  to  forecast  future 
damage  resulting  from  changes  in  the  mission  planning.  Probabilistic  modelling  using 
Bayesian  updating  has  also  been  used  to  quantify  the  potential  benefits  of  fusing 
information  from  on-board  health  monitoring  of  turbine  engine  discs.  This  analysis  fuses 
data  from  probabilistic  FaNG  model  predictions  and  continual  input  from  a  crack  detection 
sensor  to  forecast  current  and  future  probability  of  failure. 

To  facilitate  technology  transfer  to  the  U.S.  Air  Force  and  turbine  engine  community,  the 
above-mentioned  developments  that  were  judged  to  be  sufficiently  mature  and  adequately 
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validated  were  added  to  the  commercial  computer  code  Design  Assessment  of  Reliability 
With  INspection  (DARWIN®'  which  is  available  royalty-free  to  the  government  and 
through  commercial  licensing  to  industry;  these  methods,  which  are  accessible  through  the 
user  friendly  Graphical  User  Interface  (GUI),  include:  FaNG  model;  FaNG  residual  stress 
analysis  for  surface  cracks;  enhanced  mission  mixing,  and  probabilistic  sensitivities 
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2. 


Introduction 


2.1  Background 

The  USAF  is  facing  a  significant  wave  of  turbine  engine  disc  replacement  costs  in  the  next 
eight  to  ten  years  that  are  inconsistent  with  anticipated  budgets.  Consequently,  the  ERLE 
program  was  conceived  by  the  AFRL  as  a  sound  science  and  technology  investment  that 
offers  the  potential  for  significant  cost-avoidance  by  extending  the  life  of  certain  life- 
limiting  components.  The  concept  is  to  extend  the  life  of  these  components  by  recovering 
the  conservatism  believed  to  exist  in  current  design  and  life  management  practices,  without 
increasing  risk,  by  systematically  improving  and  more  effectively  integrating  a  number  of 
life  management  technologies — life  prediction,  non-destructive  inspection,  engine  health 
monitoring,  maintenance  and  repair. 

The  current  DUST  program  was  focused  on  developing  enhanced  deterministic  and 
probabilistic  methods  for  predicting  remaining  life  and  component  reliability;  the  program 
also  addressed  various  probabilistic  methods  for  enhancing  engine  health  monitoring  and 
prognosis. 

2.2  Program  Objectives 

The  overall  technical  objectives  of  the  program  were  to  develop  and  demonstrate:  1)  a 
new  family  of  physically-based,  deterministic  life  prediction  models  for  treating  total 
fatigue  life  including  crack  nucleation,  microcrack  growth  and  large  crack  growth;  2)  an 
efficient  probabilistic  life  prediction  methodology  based  on  the  differing  stochastic  nature 
of  each  of  the  above  phases  of  fatigue  life;  and  3)  an  integrated  methodology  for  enhanced 
engine  life  management  based  on  combining  and  hybridizing  state-of-the-art  probabilistic 
life  prediction  and  classical  engine  diagnostics  and  prognostics. 

The  specific  life  prediction  enhancements  pursued  in  the  current  program,  and  outlined 
below,  were  selected  from  candidate  tasks  that  were  judged  to  be  both  technically  feasible 
and  consistent  with  the  goals  of  the  ERLE  program  as  identified  in  a  prior  study  performed 
for  AFRL  by  SwRI  [1], 

2.3  Methods,  Assumptions,  and  Procedures 

The  program  was  organized  into  three  distinct  but  interrelated  technical  tasks  each  of 
which  was  aimed  at  achieving  one  of  the  above  three  program  objectives.  The  three  tasks 
were  as  follows: 

Task  1:  Enhanced  Deterministic  Life  Modelling:  This  task  developed  physically-based 
models  of  crack  nucleation,  microcrack  growth  and  large  crack  growth  phases  of  the 
fatigue  process.  These  models  were  validated  by  applying  them  to  the  prediction  of  fatigue 
damage  at  high  stress  gradients  including  those  that  can  occur  at  notches  induced  by 
foreign  object  damage  (FOD),  and  well  as  residual  stresses  resulting  from  beneficial 
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surface  treatments.  The  newly  developed  models  were  validated  by  comparing  predicted 
lives  against  available  laboratory  test  data.  Computational  simulations  were  also 
performed  to  demonstrate  the  benefits  of  tracking  fatigue  damage  based  on  actual  usage 
from  flight  data  recorders  versus  the  use  of  Total  Accumulated  Cycles  (TAC(s))  to  track 
damage  and  manage  inspection  and  maintenance  decisions.  See  Section  3.0  for  further 
details. 

Task  2:  Enhanced  Probabilistic  Life  Modelling:  This  task  utilized  the  deterministic 
models  from  Task  1  to  formulate  and  implement  a  fully-probabilistic  approach  to  total 
fatigue  life  considering  the  stochastic  nature  of  crack  nucleation,  microcrack  growth  and 
large  crack  growth.  Statistical  methods  were  also  developed  for  identifying  mission  types, 
characterizing  usage,  and  forecasting  the  effect  of  mission  mixing  changes  on  remaining 
life.  Computationally  efficient  probabilistic  methods  were  also  developed  and 
implemented  for  computing  the  effect  of  the  above  factors  on  component  reliability  so  that 
inspection,  maintenance,  and  retirement  decisions  can  be  made  on  a  reliability  basis.  See 
Section  4.0  for  further  details. 

Task  3:  Integration  of  Damage-Based  Probabilistic  Methods  with  Diagnostic  and 
Prognostics  Methods:  This  task  attempted  to  fuse  the  newly  developed  physics-based 
probabilistic  life  prediction  methods  in  Task  2  with  classical  diagnostic  and  prognostic 
methods.  The  goal  was  to  enhance  the  engine  health  management  process  by  combining 
the  rigor  and  accuracy  of  the  damage-based  life  prediction  methods  with  the  broader 
system  perspective  and  overall  efficiency  of  classical  health  monitoring  techniques  based 
on  artificial  intelligence,  neural  networks,  and  Bayesian  updating.  Bayesian  updating  was 
employed  to  demonstrate  the  benefits  of  on-board  monitoring  of  fatigue  damage  on 
component  reliability.  Limitations  were  encountered  in  attempting  to  hybridize  classical 
data  driven  methods  of  diagnostics  and  prognostics  with  physics-based  methods  and 
barriers  to  this  hybridization  are  discussed.  See  Section  4.7  for  further  details. 

2.4  Program  Team 

A  diverse  and  highly  qualified  team  was  assembled  to  address  the  technical  challenges  of 
the  proposed  program.  The  program  team  was  led  and  managed  by  Southwest  Research 
Institute,  who  served  as  prime  contractor.  Team  members  included  the  University  of 
Texas  at  San  Antonio  (UTSA),  Smiths  Aerospace,  and  Mustard  Seed  Software.  This  team 
provided  expertise  in  fundamental  understanding  of  fatigue  damage  processes  in  turbine 
engine  materials  (SwRI);  development  and  implementation  of  deterministic  and 
probabilistic  life  prediction  methods  (SwRI  and  UTSA);  ongoing  working  relationships 
with  the  major  domestic  and  foreign  engine  manufacturers  (SwRI  and  Smiths);  expertise  in 
flight  data  recording  and  engine  health  monitoring  using  on-line  data  acquisition,  data 
fusion,  and  advanced  reasoning  methods  (Smiths);  and  development  of  commercially 
successful  engineering  software  tools  for  industry  (SwRI,  Smiths,  and  Mustard  Seed 
Software). 
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2.5  Dual-Use  Benefits  and  Technology  Transfer 

The  proposed  program  directly  addresses  several  of  the  technologies  crucial  to  the  ERLE 
goal  of  doubling  the  useful  life  of  military  engines.  Studies  by  the  USAF  and  engine 
manufacturers  have  estimated  that  successful  science  and  technology  investments  of  this 
type  could  reduce  the  disk  replacement  costs  by  50%,  which  would  amount  to  a  cost 
savings  of  $600M  over  five  years.  In  addition,  enhancements  in  engine  life  management 
technology  would  also  be  applicable  to  developmental  and  future  engines.  This  is 
particularly  true  of  the  hybridization  of  life  prediction  and  health  monitoring  since  new 
engines  tend  to  contain  more  sensors  with  advanced  sensing  capability. 

Similar  benefits  would  be  expected  to  accrue  in  the  commercial  sector,  where  safety, 
reliability,  and  cost  of  ownership  are  of  paramount  importance  for  the  turbine  engines  that 
power  commercial  airliners.  In  fact,  safety  of  flight  was  the  primary  driving  force  for 
Federal  Aviation  Administration  (FAA)  and  industry  support  of  the  Turbine  Rotor  Material 
Design  (TRMD)  program,  which  led  to  the  initial  development  and  commercialization  of 
DARWIN.  The  technologies  and  associated  software  developed  in  the  current  program  are 
also  expected  to  penetrate  the  power  generation  industry  where  availability  and  cost  of 
ownership  are  high  priorities,  particularly  in  the  deregulated  market. 

Technology  demonstration  and  transfer  of  the  enhanced  probabilistic  life  prediction 
methodologies  developed  in  the  current  program  were  achieved  by  integrating  selected 
new  methodologies  as  capabilities  in  SwRI’s  DARWIN  code,  which  is  available  to 
government  agencies  under  a  royalty- free  license,  as  well  as  to  turbine  engine  and  power 
generation  communities  under  commercial  licensing. 
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3.  Results  and  Discussion  on  Enhanced  Probabilistic  Life  Modeling 

The  primary  objective  of  this  task  was  to  develop  and  demonstrate  a  new  family  of 
advanced  deterministic  life  prediction  methods  for  total  fatigue.  These  methods  would  be 
guided  by  physical  mechanisms  and  insights,  and  they  would  be  capable  of  treating  the 
various  phases  of  fatigue  life  including  crack  nucleation,  small-crack  growth,  and  large- 
crack  growth.  They  would  also  be  amenable  to  practical  engineering  use  and  statistical 
treatment.  Once  these  models  were  developed,  they  would  be  implemented  in  a  suitable 
computational  platform,  calibrated  with  available  baseline  test  data,  and  then  demonstrated 
by  comparison  with  available  feature  test  data.  These  goals  have  all  been  met,  as  described 
in  the  pages  that  follow. 

3.1  Development  of  Fatigue  Nucleation  and  Growth  (FaNG)  Model 

The  traditional  paradigm  for  fatigue  life  prediction  modeling  applied  to  engine  rotors  is 
summarized  in  Figure  1.  So-called  “crack  initiation”  analysis  is  based  on  strain-life  or 
stress-life  methods  and  is  often  calibrated  by  total  fatigue  life  results  on  small  cylindrical 
specimens.  The  end  of  the  crack  initiation  life  phase  is  commonly  assumed  to  correspond 
to  the  formation  of  an  0.030”  surface  crack.  In  order  to  account  for  normal  material 
scatter,  lower  bound  design  curves  (perhaps  corresponding  to  0.001  failure  rates)  are 
constructed  through  statistical  analysis  of  the  available  smooth  specimen  data. 

In  many  cases,  the  end  of  the  crack  initiation  phase  is  assumed  to  represent  the  end  of  the 
useful  life  of  the  component.  However,  in  some  situations,  there  is  a  desire  to  account  for 
at  least  some  fraction  of  the  useful  remaining  life  of  the  cracked  part  following  the 
initiation  event.  This  is  typically  accomplished  with  an  independent  fatigue  crack  growth 
life  computation  that  starts  with  a  0.015”  x  0.030”  surface  crack  and  ends  at  fracture.  The 
two  analyses  share  geometry  and  stress  history  information  and  some  material  properties, 
although  the  primary  materials  databases  that  support  the  independent  “initiation”  and 
“growth”  calculations  are  themselves  usually  largely  independent. 

In  this  work,  a  new  paradigm  has  been  defined  for  the  development  of  new  life  models. 

The  paradigm  is  the  simple  summation  of  crack  nucleation  and  crack  propagation  life 
phases,  with  the  transition  occurring  at  a  physically  meaningful  crack  nucleation  size.  The 
paradigm  can  be  applied  to  empirical  models,  micromechanical  models,  or  physically- 
informed  engineering  models.  The  new  paradigm  is  illustrated  generally  in  Figure  2.  Here 
history  information  may  include  stress,  strain,  and  temperature;  and  material  information 
may  include  both  traditional  empirical  mechanical  properties  and  microstructural 
information.  The  three  crack  growth  regimes  shown  may  reduce  to  two  regimes  or  even 
one  regime  for  certain  problems.  The  initial  crack  size  may  be  random  (e.g.,  defined  in 
terms  of  a  random  microstructure),  fixed,  or  dependent  on  the  nucleation  and/or  growth 
models.  Typically  the  initial  crack  size  will  be  much  smaller  than  the  traditional  0.030” 
flaw  often  used  at  the  start  of  a  conventional  damage  tolerance  analysis. 
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Figure  1.  Traditional  Paradigm  for  Engine  Rotor  Fatigue  Life  Prediction 
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Figure  2.  New  Paradigm  for  Total  Fatigue  Life  Analysis 
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This  new  paradigm  can  be  applied  to  conventional  empirical  life  models,  as  shown  in 
Figure  3.  Since  the  traditional  “crack  initiation”  models  (strain- life  or  stress-life)  estimate 
life  to  some  macrocrack  size  usually  associated  with  the  failure  of  a  smooth  laboratory 
specimen,  some  portion  of  this  life  is  actually  crack  growth.  It  is  possible  to  back-calculate 
the  crack  growth  fraction  of  the  “initiation”  life  for  some  specified  “little”  crack  size  using 
fracture  mechanics  models,  deriving  a  true  “nucleation”  curve.  A  total  crack  growth  life  to 
failure  can  then  be  calculated,  using  this  curve  as  a  starting  point.  The  total  fatigue  life  is 
then  the  summation  of  the  nucleation  and  growth  calculations,  which  are  joined  by  the 
nucleation  crack  size.  This  construction  also  permits  a  realistic  assessment  of  crack 
nucleation  vs.  crack  growth  contributions  to  total  life. 


Figure  3.  New  Total  Life  Paradigm  Applied  to  Traditional  Empirical  Models 


One  of  the  strengths  of  this  construction  is  that  it  is  possible  to  treat  directly  the  physical 
damage  growth  mechanisms  for  various  crack  sizes.  Another  strength  is  that  the  true 
nucleation  calculation  is  primarily  driven  by  the  stress  at  the  surface.  Both  of  these 
advantages  may  be  especially  relevant  when  the  application  problem  contains  significant 
stress  gradients,  as  will  be  shown. 

This  general  paradigm  was  demonstrated  and  calibrated  with  available  data  on  Ti-6A1-4V 
generated  in  the  large  USAF  High  Cycle  Fatigue  (HCF)  program  during  the  previous 
decade  [2],  The  Ti-6A1-4V  forgings  had  been  solution  treated  at  1710  °F  and  then  fan 
cooled  and  vacuum  annealed  at  1300  °F  for  2  hours.  The  resulting  microstructure 
contained  approximately  60%  primary  alpha  and  transformed  beta  lamellas.  The  yield  and 
ultimate  tensile  strengths  averaged  134.9  ksi  and  141.8  ksi,  respectively. 

Smooth  Specimen  Behavior.  A  large  number  of  fatigue  tests  on  0.25-in.  diameter  smooth 
cylindrical  specimens  were  performed  by  various  engine  companies  at  four  different  stress 
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ratios  (R  =  0.8,  0.5,  0.1,  and  -1)  as  part  of  the  Program  Research  and  Development 
Announcement  (PRDA)(industry)  program  [2,  3].  The  R  =  0.8  smooth  specimen  tests 
were  later  determined  to  exhibit  significant  time-dependent  bulk  deformation  effects,  and 
so  they  were  excluded  from  this  study.  A  modified  Smith-Watson-Topper  (SWT)  life 
model  had  been  fit  to  the  available  data  by  one  of  the  engine  companies.  This  regression 
gave 


■J<rme.mE  =  76UN  059  +  65.39  N  °'023  (]) 

where  £a,max  is  the  maximum  principal  strain  amplitude,  amax  is  the  maximum  normal  stress 
on  the  critical  plane  of  sa,max,  E  is  Young’s  modulus,  stress  and  modulus  have  units  of  ksi, 
and  N  is  total  cycles  to  failure.  The  experimental  data  and  SWT  fits  are  shown  in 
Figure  4. 

Large  Fatigue  Cracks.  All  available  Ti-6A1-4V  fatigue  crack  growth  (FCG)  rate  data 
generated  in  the  PRDA  HCF  program  were  also  collected  [2,  4],  After  tests  involving 
potentially  anomalous  load  shed  effects  were  eliminated,  a  total  of  thirty- four  large-crack 
tests  at  four  different  stress  ratios  remained  in  the  PRDA  data  base.  A  large  number  of 
tests  were  performed  at  R  =  0.1  and  R  =  0.8  (exhibiting  substantial  scatter  in  growth  rates), 
but  only  two  tests  were  conducted  and  reported  at  R  =  0.5  and  R  =  -1.  The  R  =  0.5  test 
data,  in  particular,  addressed  a  very  limited  range  of  A K. 

In  order  to  address  this  shortfall,  additional  FCG  rate  data  were  obtained  from  Ritchie  and 
Boyce  based  on  their  large-crack  investigations  in  the  Multidisciplinary  University 
Research  Initiative  MURI  (university)  HCF  program  [5],  They  had  generated  extensive 
data  at  R  =  0.1,  0.5,  and  0.8,  including  tests  at  both  50  Hz  and  1000  Hz  (no  frequency 
effects  were  observed).  The  MURI  FCG  rate  data  at  R  =  0.1  agreed  with  the  central 
tendencies  of  the  PRDA  data  but  exhibited  substantially  less  scatter.  The  MURI  data  at  R 
=  0.5  and  R  =  0.8  exhibited  slightly  slower  growth  rates  than  the  limited  PRDA  data.  Both 
the  PRDA  and  MURI  FCG  rate  data  have  also  been  statistically  analyzed  in  the  current 
program,  and  these  results  are  presented  and  discussed  in  Section  4.1. 

Tabular  representations  of  the  average  large-crack  data  suitable  for  implementation  in 
DARWIN  were  constructed.  These  tabular  representations  are  compared  with  the  original 
PRDA  and  MURI  experimental  data  in  Figures  5,  6,  and  7  for  R  =  0.1,  0.5,  and  0.8, 
respectively. 

Only  limited  fatigue  crack  growth  data — two  surface  crack  tests — were  available  for 
R  =  -1.  These  tests  were  conducted  at  significantly  different  applied  stresses  and  therefore 
exhibited  significantly  different  crack  growth  rates.  Test  “8375A”  was  conducted  at 
maximum  stresses  ranging  from  6  to  17  ksi,  and  gave  FCG  rates  ranging  from 
10'  in./cycle  to  10'  cycle.  Test  “70-22-a”  was  conducted  at  maximum  stresses  ranging 
from  58  to  70  ksi,  and  gave  growth  rates  ranging  from  around  10"6  in./cycle  to  nearly  10"4 
in./cycle. 
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Cycles  to  Failure 


Figure  4.  Smooth  Specimen  Fatigue  Life  Data  and  Smith-Watson-Topper 

Fits  for  Ti-6A1-4V 
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Figure  5.  Large-Crack  FCG  Rate  Data  and  Tabular  Model  for  Ti-6A1-4V  at  R  =  0.1 
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All  R  =  0.5  Large  Cracks 


Figure  6.  Large-Crack  FCG  Rate  Data  and  Tabular  Model  for  Ti-6A1-4V  at  R  =  0.5 
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Figure  7.  Large-Crack  FCG  Rate  Data  and  Tabular  Model  for  Ti-6A1-4V  at  R  =  0.8 


Arguments  based  on  fatigue  crack  closure  theory  suggest  that  stress  level  effects  on  fatigue 
crack  growth  rates,  which  are  minimal  at  higher  stress  ratios,  can  be  significant  at  lower 
stress  ratios.  For  example,  the  closure-based  equations  developed  by  Newman  [6]  indicate 
that  the  effective  stress  range  ratio,  U,  which  indicates  the  fraction  of  the  stress  intensity 
factor  range  that  is  effective  in  growing  the  crack,  changes  by  30%  to  50%  between  low 
stresses  typical  of  threshold  testing  and  high  stresses  typical  of  low-cycle  fatigue  (LCF) 
conditions  at  R  =  -1.  In  contrast,  U  changes  only  about  10%  at  R  =  0.1  between 
vanishingly  small  stress  levels  and  typical  LCF  stress  levels. 

Therefore,  stress-level  effects  were  considered  during  the  construction  of  representative 
tabular  curves  for  R  =  -1,  based  on  the  two  available  tests.  The  Newman  closure  equations 
were  used  to  estimate  crack  opening  levels  and  then  the  closure-corrected  range  of  the 
stress  intensity  factor  for  both  tests.  The  results  of  this  analysis  are  shown  in  Figure  8, 
along  with  a  tabular  representation  of  the  tabular  FCG  rate  curve  in  terms  of  A/Cff- 
Figure  9  then  shows  the  corresponding  representation  of  the  R  =  -1  data  in  terms  of  AK, 
along  with  hypothetical  growth  rate  curves  for  two  larger  values  of  amax.  Here  amax  = 

1 10  ksi  corresponds  to  the  largest  fully-reversed  stress  amplitude  that  will  not  introduce 
appreciable  cyclic  plasticity,  while  cmax  =  83  ksi  corresponds  approximately  to  an  LCF  life 
of  20,000  cycles  at  R  =  -1.  The  final  tabular  representation  of  the  R  =  -1  data  was 
constructed  at  a  stress  range  that  was  representative  of  the  available  smooth  specimen  and 
feature  test  data  for  R  =  -1.  All  of  the  tabular  curves  are  summarized  in  Figure  10. 
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All  R  =  -1  Large  Cracks 


Figure  8.  Closure-Corrected  FCG  Rate  Data  for  Ti-6A1-4V  at  R  =  -1 


R  =  -1 


Figure  9.  Predictions  of  FCG  Rate  Behavior  in  Ti-6A1-4V  at  R  =  -1 
Based  on  Two  Tests  and  the  Newman  Crack  Closure  Equations 
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Figure  10.  Summary  of  Tabular  Representations  for  Ti-6A1-4V 
Large-Crack  FCG  Rate  Data 


Small  Fatigue  Cracks.  Small-crack  FCG  rate  data  for  the  same  Ti-6A1-4V  were 
generated  at  R  =  0.1,  amax  =  88  ksi  by  Honeywell  (then  called  AlliedSignal)  in  the  PRDA 
program  [2,  7].  These  small-crack  data  are  compared  with  the  PRDA  large-crack  data  in 
Figure  1 1  (small-crack  data  in  red).  The  small-crack  data  exhibit  accelerated  growth  rates 
compared  to  the  large-crack  data  at  the  same  nominal  A K  values.  Also  note  that  the  small 
cracks  grew  at  nominal  A K  values  significantly  lower  than  the  apparent  large-crack 
threshold. 

A  simple  engineering  approach  to  address  small-crack  behavior  was  developed  based  on 
the  Kitagawa  [8]  diagram  construction,  Figure  12.  The  Kitagawa  diagram  integrates  the 
traditional  large  crack  fracture  mechanics  threshold,  AKth,  with  the  traditional  smooth 
specimen  endurance  limit,  Aae.  The  region  on  the  Kitagawa  diagram  both  above  the 
horizontal  endurance  limit  line  and  above  the  sloping  threshold  line  is  “unsafe”  in  that 
fatigue  failure  is  predicted  to  occur.  The  region  below  both  lines  is  “safe”  in  that  fatigue 
failure  is  predicted  not  to  occur.  The  intersection  of  the  two  lines,  defined  as  ao,  is  given 
by  the  equation 
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where,  for  consistency  with  the  smooth  specimen  endurance  limit  data  ,F=  1.12  (2/k)  is 
the  geometry  correction  factor  in  the  K  solution  for  the  semi-circular  surface  crack  in  the 
smooth  specimen. 


All  R  =  0.1  Large  Cracks  and  Small  Cracks 


Figure  11.  Large-Crack  (Multi-Colored  Symbols)  and  Small-Crack  (Red  Symbols) 
Fatigue  Crack  Growth  Rate  Data  for  Ti-6A1-4V  at  R  =  0.1 
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log  (a) 


Figure  12.  Schematic  Kitagawa  Diagram  Relating  Smooth  Specimen  Fatigue 
Endurance  Limit  (Aoe)  and  Large-Crack  Threshold  (A/Cth) 


The  Kitagawa  diagram  indicates  that  cracks  smaller  than  ao  must  be  able  to  grow  at 
nominal  stress  intensity  factors  that  are  less  than  the  large  crack  threshold,  since  smooth 
specimens  fail  by  the  initiation  and  growth  of  microcracks  to  failure  in  this  regime. 

El  Haddad  [9]  suggested  that  small  cracks  could  grow  below  the  large  crack  threshold 
because  the  effective  driving  force  of  the  small  cracks  was  larger  than  expected  from 
traditional  large  crack  fracture  mechanics.  He  proposed  that  the  actual  effective  driving 
force  was 


&Keq  =  F(a)AS^7i(a  +  a0)  ^ 

Here  F(a )  is  the  crack  shape  and  specimen/component  geometry  correction  factor  for  the 
component  of  interest;  ao  is  not  included  in  F(a).  Note  that  the  ao  contribution  is  negligible 
for  large  values  of  a,  but  becomes  increasingly  significant  for  smaller  values  of  a.  The  El 
Haddad  construction  can  also  be  interpreted  as  a  crack-size  dependent  threshold  rather  than 
a  crack-size  dependent  driving  force.  In  this  case,  the  crack-size  dependent  threshold  is 
given  as 


A£„»=A K, 


(4) 


The  El  Haddad  construction  can  be  represented  on  the  Kitagawa  diagram  as  a  dashed  line 
which  transitions  between  the  endurance  limit  and  large-crack  threshold  lines. 
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Although  the  El  Haddad  concept  has  been  criticized  for  its  apparent  lack  of  physical 
justification,  several  researchers  [10,  11]  have  independently  derived  similar  formulations 
based  on  detailed  micromechanical  considerations,  and  these  have  been  shown  to  be 
numerically  identical  to  the  El  Haddad  formulation  under  typical  conditions.  Chan  [10] 
related  ao  to  dislocation  pileups  and  mode  II  shear  cracks,  and  found  satisfactory  agreement 
with  small-crack  data  from  titanium  aluminides.  Tanaka  [11]  formulated  the  problem  in 
terms  of  a  crack-tip  slip  band  propagating  across  one  or  more  grain  boundaries,  and 
demonstrated  reasonable  agreement  between  his  formulation  and  a  variety  of  experimental 
data.  While  titanium  was  not  included  in  the  Tanaka  comparisons,  Brown  and  Taylor  [12] 
have  shown  good  agreement  with  the  El  Haddad  formulation  for  a  mill-annealed  Ti-6A1- 
4V. 

In  the  HCF  program,  this  small-crack  formulation  was  explored  primarily  from  the 
standpoint  of  crack-size  dependent  thresholds,  and  it  found  some  success  in  predicting 
fatigue  crack  behavior  associated  with  foreign-object  damage  and  fretting  fatigue  [13],  In 
the  current  program,  it  appears  to  be  more  convenient  to  apply  the  formulation  from  the 
standpoint  of  a  crack-size  dependent  driving  force,  Equation  (3),  since  the  focus  is  on 
predicting  fatigue  crack  growth  life  rather  than  threshold  stresses  (associated  with  crack 
arrest)  as  in  [13], 

The  small-crack  parameter  ao  was  estimated  to  be  approximately  0.002”  at  R  =  0.1  from 
Equation  (2)  and  the  available  experimental  data  on  the  endurance  limit  and  the  large-crack 
threshold  for  Ti-6A1-4V  generated  in  the  HCF  program.  A  similar  value  was  calculated  at 
R  =  0.5.  The  ability  of  the  crack-size  dependent  driving  force  AA7q  to  correlate  the  R  =  0.1 
small-crack  data  with  the  baseline  R  =  0.1  large-crack  data  is  shown  in  Figure  13.  The 
corrected  small-crack  data  at  even  the  smallest  crack  sizes  shown  here  (on  the  order  of 
0.001”,  or  slightly  less  than  twice  the  average  size  of  the  primary  alpha  grains)  agree  well 
with  the  large-crack  data. 
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Figure  13.  Large-Crack  (Gray  Symbols)  and  Small-Crack  (Red  Symbols) 
Fatigue  Crack  Growth  Rate  Data  for  Ti-6A1-4V  at  R  =  0.1,  Correlated  in 
Terms  of  the  El  Haddad  Crack-Size  Modified  Driving  Force  (AiCq) 


Crack  growth  rate  data  published  by  Hines,  Peters,  and  Lutjering  [14]  appear  to  provide 
some  further  evidence  of  both  the  small-crack  effect  and  the  stress-level  effect  for  Ti-6A1- 
4V.  Figure  14,  taken  directly  from  their  cited  paper,  compares  R  =  -1  small-crack  data  at 
three  different  stress  levels  with  large-crack  data  at  R  =  0.1.  The  small  cracks  exhibited 
accelerated  growth  rates  in  comparison  to  the  large-crack  data,  along  with  a  systematic 
layering  of  increasing  growth  rate  with  increasing  stress  amplitude.  Direct  analysis  of  the 
Hines  et  al.  data  is  not  possible  since  their  bi-modal  microstructure  is  slightly  different 
from  the  standard  HCF  program  material  condition,  and  since  Hines  et  al.  published  only 
idealized  trend  lines,  and  not  the  actual  data  points.  However,  it  is  possible  to  emulate  the 
Hines  et  al.  results  using  the  available  crack  growth  rate  properties  for  the  HCF  program 
Ti-6A1-4V,  and  employing  both  the  El  Haddad  small-crack  correction  as  well  as  the 
Newman  stress-level  correction.  The  resulting  “predictions”  of  the  Hines  et  al.  behavior 
are  shown  in  Figure  15  and  exhibit  the  same  general  trends.  Note  that  SI  units  are  used  in 
these  two  figures,  whereas  US  units  were  used  in  the  previous  figures. 
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Figure  14.  Hines  et  al.  (1999)  FCG  Rate  Data  for  Small  Cracks  in  Ti-6A1-4V 
Under  R  =  -1  Loading  at  Three  Different  Stress  Amplitudes,  in 
Comparison  to  Large-Crack  R  =  0.1  Behavior 


Figure  15.  Predicted  FCG  Rates  in  Ti-6A1-4V  for  R  =  -1  Conditions 
Considered  by  Hines  et  al.  (1999),  Demonstrating  Apparent 
Small-Crack  Effects  and  Stress-Level  Effects 
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FaNG  Results.  Crack  nucleation  curves  were  constructed  for  smooth  fatigue  specimens 
according  to  the  following  procedure.  First,  the  FCG  life  from  a  specified  nucleation  crack 
size  to  specimen  fracture  was  calculated  for  a  specified  uniform  stress  using  the  tabular 
large-crack  FCG  rate  data  and  the  El  Haddad  small-crack  correction  to  the  crack  driving 
force.  This  calculation  was  performed  using  the  stress  intensity  factor  for  a  surface  crack 
in  a  cylinder  (crack  case  SC07)  available  in  the  National  Aeronautics  and  Space 
Administration  (NASA)  Flaw  Grow  (NASGRO)  computer  code  [15].  Then,  this 
calculated  crack  growth  life  was  subtracted  from  the  calculated  total  life  determined  from 
the  modified  SWT  fit  to  the  PRDA  smooth  specimen  data.  The  result  was  denoted  as  the 
“nucleation  life”  associated  with  the  specified  nucleation  crack  size.  All  of  these 
algorithms  and  materials  data  were  implemented  in  a  custom  version  of  the  DARWIN 
computer  code  [16]. 

The  results  are  shown  in  Figure  16  for  R  =  0.1  and  Figure  17  for  R  =  0.5,  each  for  various 
values  of  the  nucleation  crack  size.  The  smallest  nucleation  crack  size  considered  here  is 
0.001”,  which  (as  noted  earlier)  is  slightly  less  than  twice  the  average  size  of  the  primary 
alpha  grains.  This  may  be  the  lower  limit  of  applicability  of  conventional  fracture 
mechanics  analyses.  Some  rules-of-thumb  suggest  a  slightly  larger  lower  limit,  around 
five  times  the  average  grain  diameter  (in  this  case,  about  0.003”).  The  nucleation  lives 
approach  100  percent  of  the  total  life  at  the  lowest  maximum  stresses  considered 
(approximately  the  endurance  limit)  for  both  stress  ratios.  The  nucleation  lives  are  a 
smaller  fraction  of  the  total  life  at  higher  maximum  stresses.  Note  that  the  largest 
maximum  stresses  represented  on  these  graphs  approach  the  monotonic  yield  strength  of 
the  material  (135  ksi).  The  same  results  are  represented  in  terms  of  the  relative  proportion 
of  nucleation  life  to  total  life  for  R  =  0.1  in  Figure  18.  Nucleation  life  fractions  are 
compared  at  R  =  0.1  and  R  =  0.5  for  a  fixed  nucleation  crack  size  in  Figure  19. 
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Cycles 


Figure  16.  Estimated  Crack  Nucleation  Life  and  Total  Life  for 
Ti-6A1-4V  Smooth  Specimens  at  R  =  0.1 


Figure  17.  Estimated  Crack  Nucleation  Life  and  Total  Life  for 
Ti-6A1-4V  Smooth  Specimens  at  R  =  0.5 
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Life 


Figure  18.  Estimated  Nucleation  Life  Fraction  as  a  Function  of  Total  Fatigue  Life 
at  R  =  0.1  for  Different  Nucleation  Crack  Sizes,  Calculated  Using  the  FaNG  Model 


Life 


Figure  19.  Estimated  Nucleation  Life  Fraction  as  a  Function  of  Total  Fatigue  Life 
at  R  =  0.1  and  0.5  for  a  Nucleation  Crack  Size  of  0.003  in. 
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3.2 


Model  Validation 


3.2.1  Fatigue  Crack  Growth  at  Notches 

Prediction  of  fatigue  life  at  stress  concentrations  is  a  classical  problem.  Traditional  stress- 
life  approaches  in  which  smooth  specimen  behavior  is  modified  by  the  theoretical  stress 
concentration  factor  ( kt )  give  overly  conservative  results.  Semi-empirical  corrections  to 
this  approach  employing  a  diminished  “fatigue  notch  factor”  (kf)  give  improved  results  but 
do  not  account  for  some  stress  level  and  geometry  effects. 

Better  methods  directly  address  the  partitioning  of  life  into  crack  “initiation”  and 
“propagation”  phases,  independently  calculating  each  phase  with  appropriate  methods  and 
then  summing  the  two  to  obtain  total  life.  Socie  et  al.  [17]  and  Dowling  [18]  are 
representative  notable  early  treatments  of  this  type.  The  transition  length  at  which 
initiation  and  propagation  methods  are  joined  in  these  methods  can  vary  significantly  with 
notch  geometry  and  stress  level. 

As  noted  earlier,  it  is  common  in  engine  rotor  life  management  to  associate  crack 
“initiation  life”  with  the  formation  of  a  crack  with  some  significant  (fixed)  engineering  size 
(e.g.,  0.030  in.  or  1  mm).  This  “initiation  life”  is  often  related  to  the  fatigue  failure  of  a 
smooth  cylindrical  specimen  under  uniform  axial  loading.  A  corresponding  crack 
“propagation  life”  is  then  calculated  using  this  same  size  as  the  initial  crack  size. 

However,  this  scheme  can  also  give  over-conservative  life  results,  especially  for  small 
stress  concentrations  and  steep  stress  gradients.  This  method  neglects  two  important 
phenomena:  much  of  the  “initiation  life”  is  actually  consumed  by  the  growth  of  a  very 
small  crack,  and  the  stresses  may  decrease  substantially  over  the  initiation  length. 

The  FaNG  model  demonstrated  earlier  for  smooth  specimens  can  easily  be  extended  to 
cracks  growing  at  notches  or,  more  generally,  in  stress  gradient  fields.  The  crack 
nucleation  calculation  for  the  application  problem  with  a  significant  stress  gradient  (such 
as  a  sharp  notch)  is  performed  by  determining  the  local  stresses  at  the  nucleation  length  in 
the  corresponding  (uncracked)  application  geometry,  and  then  determining  the  nucleation 
life  (at  this  nucleation  length)  for  a  smooth  specimen  subjected  to  the  same  local  stresses. 
The  local  stress  calculation  in  the  application  geometry  should  address  any  elastic -plastic 
behavior  that  may  occur  if  local  stresses  exceed  yield  at  the  stress  concentration.  The 
crack  growth  life  in  the  application  geometry  is  performed  using  customary  fracture 
mechanics  analysis  methods,  employing  stress  intensity  factor  solutions  that  address  the 
local  geometry  and  stress  gradients  in  the  application  problem.  The  same  small-crack 
correction  used  to  derive  the  nucleation  curves  must  be  used  here  to  compute  crack  growth 
life.  Finally,  the  total  fatigue  life  for  the  application  geometry  is  determined  by  summing 
nucleation  and  growth  contributions.  Because  both  life  phases  are  physically-meaningful 
values,  it  is  also  possible  to  make  a  realistic  assessment  of  the  relative  importance  of 
nucleation  and  growth  for  total  life. 
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Note  that  the  fracture  mechanics  calculation  should  not  be  performed  and  small-crack 
corrections  should  not  be  used  at  crack  sizes  so  small  that  severe  microstructural  influences 
come  into  play.  A  traditional  rule-of-thumb  is  at  least  two  times  the  grain  size,  and 
preferably  about  five  times  the  grain  size.  This  serves  to  define  a  practical  lower  limit  on 
the  nucleation  crack  size.  In  principle,  the  nucleation  crack  size  could  be  decreased  further 
if  suitable  fracture  mechanics  methods  were  available  to  characterize  the  growth  of 
microstructurally-small  fatigue  cracks.  On  the  other  hand,  it  is  not  clear  that  macroscopic 
calculations  of  the  local  stresses  (employed  in  the  nucleation  calculation)  would  be 
meaningful  at  this  microscopic  level,  and  so  it  may  be  prudent  to  maintain  the  same 
practical  lower  limit. 

The  FaNG  model  was  applied  to  predict  the  total  fatigue  lives  of  a  series  of  tests  performed 
on  double-edge-notch  (DEN)  specimens  in  the  HCF  PRD  A  program  [2],  Two  different 
notch  geometries  were  used.  The  notch  in  one  specimen  (the  so-called  “V-notch”  or 
“small  volume”  configuration)  had  a  depth  of  0.047”,  a  root  radius  of  0.021”  and  a  width 
of  0. 142”,  with  an  original  specimen  width  (neglecting  the  notch)  of  0.236”.  The  notch  in 
the  other  specimen  (“U-notch”  or  “large  volume”)  had  a  depth  and  radius  of  0.032”  and  a 
width  of  0.250”,  with  an  original  specimen  width  (neglecting  the  notch)  of  0.36”.  The 
nominal  (net-section)  stress  concentration  factor  was  approximately  2.5  for  both 
geometries.  Tests  were  performed  at  four  different  stress  ratios:  R  =  0.1,  0.5,  0.8,  and  -1. 

In  preparation  for  performing  the  FaNG  model  predictions,  three-dimensional  elastic  finite 
element  analyses  were  performed  to  obtain  the  stress  gradients  ahead  of  the  notch  root. 

The  nominal  (net-section)  stress  concentration  factor  was  approximately  2.5.  The 
maximum  elastic  stress  decreased  from  its  maximum  value  at  the  notch  root  by  a  factor  of 
two  within  0.013”  (0.33  mm)  of  the  notch  root.  The  calculated  gradients  are  shown  in 
Figure  20. 

If  the  elastic  maximum  stress  in  the  notched  specimen  was  calculated  to  be  greater  than 
yield,  then  the  DARWIN  shakedown  algorithm  was  invoked  to  calculate  the  elastic-plastic 
stress  relaxation  and  redistribution  near  the  notch  root.  When  shakedown  occurred,  the 
local  stress  ratio  was  (in  general)  different  from  the  nominal  (applied)  stress  ratio,  and  the 
actual  local  stress  ratio  was  considered  in  both  the  crack  nucleation  and  growth 
calculations.  The  cyclic  stress-strain  properties  of  Ti-6A1-4V  were  employed  in  the 
shakedown  calculations;  the  cyclic  yield  strength  is  about  1 13  ksi.  An  example  of  the 
original  elastic  and  shakedown  elastic-plastic  stress  gradients  is  provided  in  Figure  21  for 
nominal  R  =  0.5  loading  with  a  maximum  nominal  applied  stress  of  57.5  ksi.  This  figure 
also  illustrates  the  definition  of  the  “nucleation  stress”  at  the  specified  nucleation  length  (in 
this  case,  0.003”). 
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Figure  20.  Calculated  Elastic  Stress  Gradients  for 
Two  Double-Edge-Notch  Specimen  Geometries 
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Figure  21.  Calculated  Elastic  and  Elastic-Plastic  Shakedown  Stress  Gradients 
at  Maximum  and  Minimum  Load  for  Nominal  R  =  0.5  Loading 
with  a  Remote  Applied  Stress  of  57.5  ksi 
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The  test  reports  indicated  that  most  of  the  observed  cracking  initiated  as  a  dominant  semi¬ 
elliptical  surface  crack  near  the  center  of  the  notch  root,  so  the  initial  conditions  for  the 
crack  growth  analysis  were  selected  as  a  semi-circular  crack  (with  depth  equal  to  the 
nucleation  length)  at  mid-section.  The  DARWIN  stress  intensity  factor  solution  SC  17,  a 
weight  function  solution  for  a  semi-elliptical  surface  crack  growing  in  a  univariant  stress 
gradient  field,  was  used  for  the  calculations. 

Several  different  values  of  the  nucleation  length  were  considered  during  early 
investigations.  Note  that  as  the  nucleation  length  is  decreased,  the  nucleation  life 
associated  with  that  length  decreases,  but  the  crack  growth  life  to  failure  associated  with 
that  nucleation  length  increases.  The  calculated  total  life  (nucleation  +  growth)  will  not 
change  with  nucleation  length  in  the  absence  of  a  stress  gradient.  When  the  local 
nucleation  stress  changes  significantly  with  nucleation  length,  however,  the  calculated  total 
life  will  change  with  nucleation  length. 

Comparisons  of  model  predictions  for  total  life  with  DEN  test  data  at  R  =  0.1  for  three 
different  values  of  the  nucleation  length  are  shown  in  Figure  22.  The  three  different  values 
of  nucleation  length  give  similar  total  life  predictions  at  higher  maximum  stresses,  because 
yielding  near  the  notch  root  (which  causes  stress  relaxation  and  redistribution) 
substantially  reduces  the  severity  of  the  stress  gradient.  At  lower  maximum  stresses,  the 
different  nucleation  lengths  had  a  larger  effect  on  the  predicted  total  life,  although  the 
differences  in  predicted  life  were  still  sometimes  within  the  scatter  of  the  experimental  life 
results.  A  nucleation  length  of  about  0.003”  generally  gave  the  most  accurate  predictions. 

Also  shown  in  Figure  22  for  reference  purposes  are  total  life  predictions  based  on  a 
conventional  stress-life  approach  in  which  the  modified  SWT  model  is  applied  to  the 
elastically-calculated  stresses  and  strains  at  the  notch  root  itself  (essentially,  the  SWT 
model  is  scaled  by  the  peak  stress  concentration  factor).  As  noted  previously,  this  type  of 
simple  approach  substantially  overestimates  the  fatigue  damage  caused  by  the  notch  and 
gives  over-conservative  total  life  estimates.  Plastic  corrections  to  the  local  stresses  and 
strains  (non-negligible  only  above  45  ksi)  would  be  inadequate  to  reconcile  the  differences. 

Figure  23  shows  comparisons  of  predictions  and  test  data  at  R  =  0.1  for  the  two  notch 
geometries,  considering  only  the  0.003”  nucleation  length.  The  two  geometries  are  shown 
to  be  quite  similar,  exhibiting  only  minor  life  differences  in  both  test  and  model. 

Model  predictions  and  test  data  for  R  =  0.5  (all  employing  the  V-notch  geometry)  are 
compared  in  Figure  24.  Again,  the  normalized  SWT  model  is  very  conservative.  Similar 
results  are  shown  for  R  =  -1  conditions  (all  V-notch)  in  Figure  25. 


26 


R  =  0.1  (V-Notch) 


l.E+04  l.E+05  l.E+06  l.E+07  l.E+08 

Life  (cycles) 


Figure  22.  Comparison  of  FaNG  Model  Total  Life  Predictions  and  Test  Data 
for  V-Notch  Specimens  at  R  =  0.1  with  Different  Nucleation  Sizes 
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Figure  23.  Comparison  of  FaNG  Model  Total  Life  Predictions  and  Test  Data 
for  V-Notch  and  U-Notch  Specimens  at  R  =  0.1 
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Figure  24.  Comparison  of  FaNG  Model  Total  Life  Predictions  and  Test  Data  at  R  =  0.5 
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Figure  25.  Comparison  of  FaNG  Model  Total  Life  Predictions  and  Test  Data  at  R  =  -1 
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FaNG  model  predictions  of  the  V-notch  fatigue  data  for  all  three  stress  ratios  are 
summarized  in  Figure  26,  considering  only  the  0.003”  nucleation  length.  Agreement 
between  test  and  model  is  generally  strong  under  all  conditions.  It  should  be  emphasized 
that  the  predictions  shown  here  are  in  no  way  a  fit  to  the  notch  test  data  (other  than  the 
selection  of  a  common  nucleation  length),  because  the  material  constants  are  derived  only 
from  smooth  specimen  and  fracture  mechanics  tests. 

The  model  facilitates  evaluation  of  the  fraction  of  total  fatigue  life  consumed  by  crack 
nucleation,  for  the  specified  nucleation  length.  These  predictions  are  summarized  in 
Figure  27  for  the  same  three  stress  ratios  and  the  same  nucleation  length  (0.003”).  In  the 
long  life  regime,  the  total  life  is  dominated  by  crack  nucleation,  but  in  the  low  cycle  regime 
(total  lifetimes  of  105  cycles  and  less),  the  crack  growth  phase  consumes  30%  and  more  of 
the  total  life.  It  is  interesting  to  note  that,  for  this  particular  geometry,  the  crack  nucleation 
life  fraction  for  a  given  total  life  does  not  change  significantly  with  stress  ratio. 

This  study  must  be  regarded  as  a  preliminary  investigation  due  to  its  limited  scope. 

Further  work  is  needed  to  evaluate  the  method  for  other  materials  and  other  feature 
geometries.  Selection  of  the  optimum  “nucleation  length”  also  requires  further  study.  It  is 
not  yet  clear  if  this  length  is  dependent  on  material  parameters,  geometry  parameters,  stress 
gradient  parameters,  or  some  combination  of  the  three.  Although  the  nucleation  length 
does  not  always  appear  to  play  a  significant  role  in  the  quality  of  the  life  prediction,  it  is 
not  yet  possible  to  make  an  a  priori  determination  of  the  only  free  parameter  in  the 
analysis  scheme. 

Nevertheless,  the  results  of  this  first  study  are  encouraging.  In  principle,  the  method 
should  apply  not  only  to  stress  concentrations,  but  also  to  other  problems  with  significant 
stress  gradients.  For  example,  many  surface  enhancement  methods  (such  as  peening) 
induce  a  severe  stress  gradient  very  close  to  the  surface.  These  applications  will  be 
explored  in  the  following  section.  The  new  method  provides  a  rational  basis  to  partition 
the  problem  into  physically-meaningful  nucleation  and  growth  phases,  to  assess  the 
relative  importance  of  each,  and  to  obtain  improved  life  predictions. 
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Figure  26.  Summary  Comparison  of  FaNG  Model  Total  Life  Predictions 
and  Test  Data  at  Three  Stress  Ratios  for  0.003”  Nucleation  Size 
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Figure  27.  Predicted  Nucleation  Life  Fraction  for  Notched  Geometry 
and  0.003”  Nucleation  Size  at  Three  Stress  Ratios 


30 


3.2.2  Near-Surface  Residual  Stresses 


Another  significant  class  of  fatigue  problems  involving  significant  stress  gradients  at  the 
surface  is  the  general  category  of  surface  enhancements  such  as  shot  peening,  laser  shock 
peening,  or  low  plasticity  burnishing.  These  surface  enhancements  are  applied  for  the 
specific  purpose  of  inducing  a  pronounced  compressive  residual  stress  at  the  surface.  This 
residual  stress,  along  with  the  associated  local  plastic  deformation,  generally  leads  to  a 
substantial  increase  in  fatigue  life.  However,  analysis  of  these  benefits  is  difficult  because 
of  the  severity  of  the  stress  gradients  over  very  small  length  scales. 

The  FaNG  model  appears  to  be  an  appropriate  analytical  choice  to  address  this  class  of 
problems.  However,  the  severity  of  the  local  stress  field  introduces  some  practical 
difficulties  in  performing  the  fatigue  crack  growth  calculation  using  weight  functions.  The 
local  stress  values  in  these  problems  can  change  dramatically  and  even  reverse  sign  over 
very  small  distances.  At  the  same  time,  the  component  may  contain  an  applied  stress  field 
with  a  gradient  operative  over  a  much  longer  length  scale.  These  stress  fields  with 
dissimilar  length  scales  present  a  computational  challenge  to  the  life  prediction  calculation. 
Highly  irregular  point  spacing  can  compromise  computational  efficiency,  which  is 
especially  critical  for  probabilistic  analysis.  For  this  reason,  the  DARWIN  weight  function 
solutions  were  originally  configured  for  stresses  defined  at  equally-spaced  points  along  the 
gradient  line.  This  approach  was  feasible  for  applications  problems  in  which  stress 
gradients  exhibited  a  single  length  scale,  but  the  approach  breaks  down  for  problems 
involving  multiple  gradients  with  multiple  length  scales.  Sparse  point  spacing  appropriate 
for  a  large  length  scale  does  not  provide  enough  resolution  to  accurately  define  a  local 
gradient  over  a  small  length  scale;  but  fine  spacing  appropriate  for  the  local  gradient  will 
require  an  excessively  large  number  of  points  for  the  gradient  over  the  large  length  scale. 

In  order  to  address  this  problem,  a  new  optimum  point  spacing  (OPS)  algorithm  was 
developed  for  weight  function  stress  intensity  factor  solutions  in  DARWIN.  The  goal  of 
the  algorithm  is  to  define  the  smallest  number  of  points  and  the  optimum  spacing  of  those 
points  to  represent  stress  variations  with  high  fidelity  and  without  impacting  the  efficiency 
of  search  routines  in  the  stress  intensity  factor  solutions.  The  OPS  algorithm  uses  Hermite 
interpolation  to  define  a  reference  curve  which  is  optimized  such  that  the  OPS  and  user- 
defined  stresses  agree  within  0.1%.  The  point  density  of  the  OPS  representation  increases 
at  locations  with  sudden  change  in  stress  gradient. 

Figure  28  (left-hand  side)  illustrates  the  limitations  of  equal  point  spacing.  The  point 
spacing  is  insufficient  to  capture  the  actual  stress  variation  near  the  origin,  but  excessive  in 
describing  the  constant  stress  field  farther  away  from  the  origin.  The  right-hand  side  of 
Figure  28  demonstrates  how  the  OPS  algorithm  selects  29  points  to  describe  the  entire 
curve  accurately  and  efficiently. 


31 


Normalized  Distance  Normalized  Distance 


Figure  28.  Illustration  of  Equal  Point  Spacing  (Left) 
and  Optimum  Point  Spacing  Algorithm  (Right) 


Fatigue  data  [19]  provided  by  Lambda  Research  were  available  for  coupons  machined 
from  the  same  Ti-6A1-4V  material  used  in  the  HCF  program  and  then  subjected  to  shot- 
peening  (SP)  or  low-plasticity  burnishing  (LPB).  Following  surface  enhancement  but  prior 
to  fatigue  testing,  the  coupons  were  also  thermally  exposed  for  10  hours  at  795  °F 
(425  °C).  Residual  stress  measurements  were  made  on  similar  specimens  after  peening, 
both  before  and  after  thermal  exposure.  The  residual  stress  measurements  following 
thermal  exposure  are  shown  in  Figure  29.  The  thermal  exposure  caused  significant 
relaxation  of  the  original  SP  residual  stresses,  but  relatively  little  relaxation  of  the  original 
LPB  residual  stresses.  The  fatigue  tests  were  conducted  at  R  =  0.1. 

As  shown  in  Figure  30,  a  comparison  of  the  Lambda  fatigue  test  results  with  the  original 
HCF  PRDA  smooth  specimen  stress  versus  number  of  cycles  to  failure  (S-N)  data  revealed 
that  the  observed  fatigue  lives  for  the  SP  specimens  were  actually  shorter  than  the  observed 
fatigue  lives  for  the  baseline  (unpeened)  specimens  tested  in  the  industry  HCF  program 
(these  are  the  tests  that  were  used  earlier  in  this  program  to  calibrate  the  FaNG  model). 

The  LPB  specimens  exhibited  slightly  higher  fatigue  strengths  than  the  original  baseline 
specimens  from  the  industry  program.  Lambda  Research  has  not  published  any  baseline 
(unpeened)  data  of  their  own  for  this  material  under  similar  testing  conditions  (four-point 
bend  with  a  trapezoidal  cross-section).  Not  included  in  this  figure  are  runout  (usually 
around  2E6  cycles)  results  (LPB  tests  with  max  stresses  from  85  to  93  ksi,  and  SP  tests 
with  max  stresses  from  65  to  70  ksi);  the  only  data  shown  here  are  the  actual  failures. 
Lambda  specimens  with  FOD  damage  have  also  been  excluded  from  consideration  here. 
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Figure  29.  Measured  Residual  Stress  Profiles  in  Ti-6A1-4V  Coupons  for  Shot-Peening 
(SP)  and  Low-Plasticity  Burnishing  (LPB)  Followed  by  Thermal  Exposure 
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Figure  30.  Comparison  of  Lambda  Data  for  SP  and  LPB  Specimens 
with  HCF  PRDA  Baseline  Smooth  Specimen  Fatigue  Data 
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Preliminary  studies  with  the  conventional  FaNG  model  found  that  the  severe  compressive 
stresses  associated  with  the  engineered  residual  stresses  led  to  calculated  nucleation 
lifetimes  that  were  much  longer  than  the  observed  total  fatigue  lifetimes;  in  fact,  the 
calculated  nucleation  lifetimes  were  infinite  in  many  cases,  because  the  local  stresses  were 
below  the  smooth  specimen  endurance  limit. 

Therefore,  it  was  postulated  that  the  crack  nucleation  lifetime  was  zero  for  both  SP  and 
LPB.  This  might  be  consistent,  for  example,  with  heavy  shot  peening  causing  damage  to 
the  specimen  surface  and  resulting  in  immediate  crack  formation.  The  physical  motivation 
for  this  assumption  is  less  clear  for  LPB,  but  it  is  interesting  to  explore  the  implications  of 
the  assumption.  The  usual  FaNG  crack  growth  life  calculation  was  then  performed  using 
the  same  tabular  fatigue  crack  growth  data,  the  same  small-crack  model,  the  same  weight 
function  stress  intensity  factor  approach,  and  the  same  crack  nucleation  length  (initial 
crack  length)  as  in  the  previous  validation  study  on  notched  Ti-6AL-4V  specimens  from 
the  HCF  PRDA  program.  The  measured  residual  stress  field  (SP  followed  by  thermal 
exposure)  reported  by  Prevey  et  al.  was  superimposed  with  a  uniform  tensile  applied  stress 
at  R  =  0. 1 .  The  actual  Lambda  geometry  is  a  four-point  bend  specimen,  but  the  differences 
between  a  bend  or  uniform  tension  field  are  negligible  at  the  shortest  crack  sizes  that 
control  the  fatigue  life.  The  OPS  algorithm  was  used  to  superimpose  the  stress  fields. 

The  calculated  local  stress  ratios  for  the  smallest  cracks  in  the  presence  of  residual  stresses 
are  sometimes  strongly  negative  (lower  than  R  =  -1).  However,  baseline  FCG  rate  data  are 
not  available  at  stress  ratios  lower  than  R  =  - 1 .  In  order  to  perform  FCG  analyses  at  these 
low  stress  ratios,  tabular  FCG  rate  curves  at  R  =  -2,  -3,  and  -5  were  derived  from  the 
available  R  =  -1  data  using  the  Newman  equations  [6]  (based  on  crack  closure  theory)  for 
the  effect  of  stress  ratio. 

Figure  31  compares  the  Prevey  et  al.  Ti-6A1-4V  HCF  data  for  both  SP  and  LPB  conditions 
with  the  baseline  smooth  specimen  data  from  the  PRDA  program  and  the  predictions  of  the 
FaNG  crack  growth  model.  The  agreement  between  the  test  data  and  the  FCG  predictions 
is  striking.  It  should  be  emphasized  that  these  predictions  are  entirely  independent  of  the 
test  results  against  which  they  are  compared  in  this  figure;  there  is  no  fitting  or  calibration 
involved,  and  all  of  the  material  properties  used  in  the  predictions  were  derived  from 
completely  independent  tests.  The  implication  of  these  results  is  that,  at  least  for  these 
coupons,  the  crack  nucleation  life  (to  these  small  initial  crack  sizes)  appears  to  be 
negligible. 

Although  the  LPB  predictions  agreed  well  with  the  central  tendencies  of  the  test  data,  the 
test  data  exhibited  considerable  scatter  about  this  central  tendency,  and  further  studies  were 
conducted  to  investigate  the  probable  sources  of  this  scatter. 
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Figure  31.  Comparisons  of  Experimental  Fatigue  Life  Data  and  FCG  Life  Predictions 
for  SP  and  LPB  Coupons  After  Thermal  Exposure 
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In  the  first  sensitivity  study,  the  initial  crack  size  (the  nucleation  crack  size)  was  varied 
between  0.002”  and  0.006”.  The  resulting  family  of  predicted  life  curves  is  shown  in 
Figure  32.  Although  the  small  variability  in  life  for  these  different  nucleation  crack  sizes 
is  consistent  with  the  small  variability  in  observed  test  life  for  the  SP  coupons,  the 
variability  in  test  life  for  the  LPB  coupons  is  generally  larger  than  would  be  explained  by 
reasonable  variability  in  initial  crack  sizes.  Therefore,  other  sources  of  variability  were 
investigated. 

Lambda  Research  has  published  another  study  [20]  of  thermal  residual  stress  relaxation 
including  Ti-6A1-4V.  In  this  study,  the  residual  stress  profiles  were  measured  for  both  SP 
and  Laser  Shock  Peened  (LSP)  coupons  immediately  after  surface  enhancement  and  then 
after  various  times  of  thermal  exposure  at  615  °F  (325  °C).  The  residual  stress 
measurements  are  shown  in  Figure  33  for  SP  coupons  and  Figure  34  for  LSP  coupons. 
Although  there  is  clearly  a  relaxation  in  the  residual  stresses  between  the  baseline 
(unexposed)  and  thermally  exposed  coupons,  there  does  not  appear  to  be  a  consistent 
change  in  the  residual  stress  profile  for  progressively  longer  exposure  times.  Although  a 
progressive  decrease  in  the  residual  stresses  with  longer  exposure  times  is  physically 
reasonable,  it  is  possible  that  it  is  confounded  by  measurement  error  or  by  location-to- 
location  variability  in  the  initial  residual  stress  field  (no  replicate  measurements  were 
published  in  this  study).  The  logical  next  question  is  how  much  impact  this  variability  in 
the  residual  stresses  themselves  might  have  on  fatigue  life  variability  for  treated  coupons. 
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Figure  32.  Comparisons  of  Experimental  Fatigue  Life  Data  for  SP  and  LPB  Coupons 
After  Thermal  Exposure  with  Fatigue  Crack  Growth  Life  Predictions 
Using  Different  Initial  Crack  Sizes 


Figure  33.  Measured  Residual  Stresses  for  Ti-6A1-4V  Test  Coupons  After  Shot  Peening 
(Baseline)  and  After  Subsequent  Thermal  Exposure  (615  °F)  for  Various  Lengths  of  Time 
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Figure  34.  Measured  Residual  Stresses  for  Ti-6A1-4V  Test  Coupons  After 
Laser  Shock  Peening  (Baseline)  and  After  Subsequent  Thermal  Exposure 
(615  °F)  for  Various  Lengths  of  Time 


A  simple  sensitivity  study  was  conducted  to  address  this  question.  It  was  observed  that  the 
various  measured  residual  stress  profiles  in  Figures  33  and  34  could  be  enveloped  by 
“bounding”  curves  that  deviated  from  the  “central  tendency”  curve  by  either  ±10  ksi  (SP) 
or  ±8  ksi  (LPB).  Here  the  central  tendency  and  bounding  curves  were  informally 
constructed  by  eye,  not  with  any  formal  regressions.  Then  the  residual  stress 
measurements  shown  in  Figure  29  (and  used  to  generate  the  life  predictions  in  Figure  30) 
were  arbitrarily  adjusted  by  shifting  the  entire  curve  up  or  down  by  9  ksi.  Finally,  the 
adjusted  residual  stress  curves  were  used  to  perform  similar  fatigue  crack  growth  life 
predictions,  using  the  same  nucleation  crack  size  employed  earlier  (0.003”).  The  resulting 
predictions  are  shown  in  Figure  35.  The  variability  in  fatigue  life  observed  in  the 
experiments  is  similar  to  the  life  variability  predicted  by  this  relatively  minor  variability  in 
residual  stress.  This  is  no  proof,  of  course,  but  it  suggests  that  further  investigations  of 
residual  stress  variability  and  its  implications  for  fatigue  life  variability  may  be  warranted. 

Again,  this  study  must  be  regarded  as  a  preliminary  investigation  due  to  its  limited  scope. 
Further  work  is  needed  to  evaluate  the  method  for  other  materials  and  other  residual  stress 
fields.  However,  these  first  results  are  encouraging  and  suggest  that  a  careful  treatment  of 
fatigue  crack  growth,  following  well-defined  methodologies,  offers  promise  for  accurate 
evaluation  of  surface  enhancement  effects  on  fatigue  life.  Other  applications  involving 
steep  stress  gradients,  such  as  those  generated  by  large  thermal  gradients,  should  also  be 
explored. 
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Figure  35.  Comparisons  of  Experimental  Fatigue  Life  Data  for  SP  and  LPB  Coupons 
Tested  After  Thermal  Exposure  with  FCG  Life  Predictions  Using  the 
Measured  Residual  Stress  Profile  Shifted  Up  or  Down  by  9  ksi 
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3.3  Predicting  Fatigue  Life  from  TACs  Versus  Usage 

The  US  Department  of  Defense  has  traditionally  used  the  concept  of  TACs  for 
maintenance  and  inspection  scheduling,  fatigue  life  management,  and  life  warranty 
assessment  of  fatigue  and  fracture  critical  turbine  engine  components.  This  concept 
attempts  to  characterize  the  often  complex  engine  usage  by  a  number  of  effective  single¬ 
amplitude,  LCF  cycles.  The  basic  definition  of  a  TAC  is  provided  by  an  industry 
consensus  Society  of  Automotive  Engineering  (SAE)  International  standard  as  follows 
[21]: 


TAC  =  LCF  +  FTC/4  +  PTC/40  (5) 

where 

LCF  =  Low-cycle  fatigue  cycles  (“Intermediate”  and  above  and  back  to  “off’) 

FTC  =  Full-throttle  cycle  (“Intermediate”  and  above  and  back  to  “idle”) 

PTC  =  Part-throttle  cycle  (“Cruise”  to  Intermediate  and  above  and  back  to  “cruise”) 

and  the  specific  revolutions  per  minute  (RPM)  values  for  these  terms  are  engine  specific. 
As  can  be  seen  from  Equation  (5),  the  FTC  cycles  are  assumed  to  be  !4  as  damaging  as  the 
LCF  cycle,  consequently  four  of  these  FTC  cycles  would  be  required  to  produce  damage 
equivalent  to  the  single  LCF  cycle.  Analogously,  40  of  the  much  small  amplitude  PTC 
cycles  would  be  required  to  produce  damage  equivalent  to  the  single  LCF  cycle.  Equation 
(5)  was  developed  by  comparing  computed  damage  from  actual  measured  usage  with  that 
from  TACs  using  classical  fatigue  “crack  initiation”  analyses,  where  “crack  initiation”  is 
defined  as  the  failure  of  a  smooth,  approximately  0.25  in.  diameter,  fatigue  specimen.  In 
other  words,  the  best  constants  (1/4  and  1/40)  in  Equation  (5)  were  calibrated  using  the 
above  analyses  for  a  usages  from  a  variety  of  different  missions. 

The  work  summarized  in  the  remainder  of  this  section  is  based  on  a  similar  analysis  using 
several  different  fatigue  models  including  classical  “crack  initiation”,  growth  of  small  and 
large  cracks,  as  well  as  selected  combinations  of  the  above,  as  described  in  more  detail 
later.  The  object  of  this  study  was  develop  a  better  understanding  of  fatigue  lives 
computed  from  TACs,  as  well  as  to  assess  whether  or  not  lives  computed  from  TACs  are 
in  general  conservative  relative  to  those  computed  from  actual  mission  usage  experienced 
by  engines. 

3.3.1  Mission  Usages  Considered 

Actual  mission  usage  time  histories  were  provided  from  Flight  Data  Recorders  (FDRs)  by 
Smiths  Aerospace.  Since  these  FDRs  include  algorithms  that  counts  events  (i.e.,  LCF, 
FTC,  PTC)  associated  with  TACs  based  on  RPM-level  crossings,  the  algorithms  employ  a 
modified  expression  for  TACs  to  account  for  level  crossings  that  occur  within  multiple 
events.  For  example,  an  LCF  cycle  includes  also  the  level  crossings  associated  with  a  FTC 
and  a  PTC,  but  should  be  counted  as  a  single  TAC.  The  modified  TAC  expression  that 
takes  this  into  account  is  as  follows: 
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TAC  =  LCF  +  (FTC  -  LCF)/4  +  (PTC  -  FTC)/40 


(6) 


where  for  the  engine  of  interest  here  the  events  are  defined  in  terms  of  either  low  speed 
(Nl)  or  high  speed  (N2)  spools  of  the  engine  as  follows: 

N1  Spool: 

LCF  -  RPM  below  2000  to  above  9500  and  back  below  2000 
FTC  -  RPM  below  5900  to  above  9500  and  back  below  5900 
PTC  -  RPM  below  7700  to  above  9500  and  back  below  7700 

N2  Spool: 

LCF  -  RPM  below  6500  to  above  12500  and  back  below  6500 
FTC  -  RPM  below  10250  to  above  12500  and  back  below  10250 
PTC  -  RPM  below  11500  to  above  12500  and  back  below  11500 

An  application  example  problem  was  developed  to  investigate  the  relationship  between 
TACs  and  the  number  of  flights  to  failure.  Stress  values  were  obtained  by  applying  N2 
(high  speed)  RPM  values  from  FI 00  usage  histories  provided  by  Smiths  Aerospace  to  a  3D 
finite  element  engine  disk  model  calibrated  to  reference  stress  and  RPM  values.  Nearly 
200  usage  histories  were  considered  initially.  However,  many  of  the  histories  were 
eliminated  for  various  reasons  (e.g.,  mission  type  could  not  be  identified,  usage  history  was 
incomplete  or  was  not  an  actual  airborne  mission).  Of  the  remaining  3 1  histories, 
indicated  in  Table  1,  most  of  the  data  were  in  the  Live  Fire  and  Instruments  &  Navigation 
categories. 

Table  1.  F100  Engine  Usages  Considered  for  Crack  Growth  Life  Comparisons 


Mission  Type 

Elmendorf 

Lakenheath 

Mtn.  Home 

Total 

Live  Fire 

4 

0 

4 

8 

Live  Fire  Missiles 

0 

2 

1 

3 

ACBT 

0 

2 

0 

2 

Air  to  Ground  Weapons 

0 

1 

1 

2 

Instruments  &  Navigation 

0 

8 

5 

13 

Transition  &  Proficiency 

0 

0 

2 

2 

Functional  Check  Flight 

0 

1 

0 

1 

Total 

4 

14 

13 
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3.3.2  Computed  Fatigue  Lives  from  Actual  Mission  Usage  Versus  TACs 

Crack  growth  life  and  TAC  values  were  computed  for  a  representative  disk  using  stresses 
based  on  the  N2  RPM  values  from  FI 00  usage  histories  provided  by  Smiths  Aerospace 
(See  Table  1). 
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To  facilitate  comparison  of  computed  results  for  different  mission  usages  a  normalization 
constant  Nf*  was  selected  based  on  the  crack  growth  life  associated  with  the  maximum  (0- 
100%)  design  RPM  for  the  FI 00  engine  (14,300  RPM  for  N2  spool),  combined  with  the 
number  of  TACs  per  flight,  as  shown  below: 

*  _  Flights  to  Failure  at  100%  Design  RPM  (7) 

f  TACs  per  flight 

Fatigue  life  predictions  were  performed  using  several  different  fatigue  models  as  follows: 

1)  Classical  “crack  initiation”  using  the  SWT  parameter 

Nf  =  N;SWT  (8) 


2)  Classical  long-crack  fracture  mechanics 


Nf  =  NLC  (9) 

3)  SWT  “crack  initiation”  plus  larger  crack  growth  from  a  0.030  in.  “initiation”  crack  size 

Nf=NiSWT  +  NLc  (10) 

4)  The  “true”  FAtigue  Nucleation  and  Growth  model  described  in  Section  3.1 


Nf  -  Ntn  +  (Nsc  +  Nlc) 


(11) 


where  N  represents  the  number  of  cycles  and  superscript  SWT  =  SWT  parameter  and 
subscripts  f  =  failure,  i  =  initiation,  tn  =  true  nucleation,  LC  =  large  crack,  and  SC  =  small 
crack. 

Figure  36  shows  results  based  on  classical  crack  initiation  and  the  SWT  parameter.  As  can 
be  seen  all  of  these  normalized  results  are  greater  than  1  indicating  the  use  of  TACs  is 
conservative  by  at  least  a  factor  of  two  on  life.  This  is  likely  due  to  the  fact  that  the  mean 
life  curve  of  SWT  versus  Nf  was  employed  in  these  calculations,  whereas  using  a  mean 
minus  standard  deviation  curve  would  reduce  to  lives. 

In  comparison  Figure  37  shows  results  based  on  the  classical  long-crack  fracture 
mechanics  (Equation  (8))  with  an  assumed  initial  surface  flaw  size  of  15  x  30  mils.  As  can 
be  seen,  many  of  the  normalized  life  values  fall  below  one  indicating  that  results  can  be 
non-conservative  when  using  TACs  to  compute  fatigue  life  for  this  case.  However,  as  one 
might  expect  the  behavior  changes  significantly  when  nucleation  plus  long  crack  growth  is 
considered  as  shown  in  Figure  38;  these  results  are  equivalent  to  adding  the  results  from 
Figures  36  and  37.  Now  the  predicted  lives  based  on  TACs  are  again  conservative  and  in 
fact  are  very  similar  to  the  results  from  the  initiation  only  calculations  shown  in  Figure  36. 
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This  similarity  is  due  to  the  fact  that  when  one  assumes  a  relatively  long  transition  crack 
size  between  initiation  and  growth  of  15  X  30  mils,  most  of  the  fatigue  life  is  dominated  by 
“initiation”.  This  is  the  reason  for  the  commonly  held  perception  in  the  turbine  engine 
industry  that  LCF  of  discs  is  an  “initiation”  event. 

Figure  39  compares  the  SWT  initiation  lives  with  the  true  nucleation  lives  as  defined 
within  the  context  of  the  FaNG  model.  As  can  be  seen,  the  resulting  livers  from  the  two 
models  are  noticeably  different  since  the  true  nucleation  lives  are  for  a  nucleation  crack 
depth  of  3  mils  (see  Section  3.2)  versus  that  corresponding  to  failure  of  a  smooth  S-N 
specimen,  which  is  closer  to  30  mils. 

If  we  now  compare  the  total  fatigue  lives  computed  from  Equations  (10)  versus  (1 1)  as 
shown  in  Figure  40  we  see  that  results  are  very  similar.  This  is  because  the  stress  gradient 
in  this  problem  is  rather  mild  and  thus  the  definition  of  the  transition  crack  size  does  not 
measurably  affect  the  computed  fatigue  lives.  However  if  we  tried  to  compare  results  from 
these  same  two  models  for  a  strip  stress  gradient,  the  results  would  be  significantly 
different  since  the  classical  initiation  crack  size  of  30  mils  would  likely  be  well  outside  the 
region  of  stress  concentration.  This  is  precisely  why  the  classical  fatigue  initiation 
methods  have  difficulty  with  correctly  predicting  notch  effects  particularly  when  the 
gradient  (or  notch  root  radius)  varies  from  one  case  to  the  next. 

It  is  interesting  to  note  that  for  the  missions  considered,  as  well  as  for  most  of  the  models, 
the  normalized  life  values  are  significantly  greater  than  one  and  increase  with  an  increasing 
number  of  TACs.  This  suggests  that  the  TACs-based  relationship  is  conservative  when 
applied  to  total  life  values,  and  conservatism  increases  as  the  number  of  TACs  increases. 
Consequently,  the  benefit  from  employing  actual  usage,  rather  than  TACs,  in  tracking 
remaining  fatigue  life  is  that  some  of  this  conservative  margin  in  the  calculations  could  be 
used  to  extend  the  component  life. 
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Figure  36.  Normalized  Nucleation  Life  (Smith-Watson-Topper  Model)  Associated 
with  F100  Engine  Usage  Versus  Total  Accumulated  Cycles  per  Flight 
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Figure  37.  Normalized  Large  Crack  Propagation  Life  Values  Associated  with  F100 
Engine  Usage  Versus  Total  Accumulated  Cycles  per  Flight 
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Figure  38.  Normalized  Nucleation  (Smith- Watson-Topper)  Plus  Propagation 
Life  Associated  with  F100  Engine  Usage  Versus 
Total  Accumulated  Cycles  per  Flight 


Figure  39.  Normalized  Crack  Nucleation  Lives  Associated  with 
Smith-Watson-Topper  and  FaNG  Models 
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Figure  40.  Normalized  Crack  Nucleation  Plus  Propagation  Lives  Associated 
with  Smith-Watson-Topper  and  FaNG  Models 
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4.  Results  and  Discussion  on  Enhanced  Probabilistic  Life  Modeling 

4.1  Modeling  Uncertainty  in  Fatigue  Crack  Growth  Rates 

The  dependence  of  fatigue  crack  growth  rates  on  AK  and  R-value,  da/dN(AK,  R),  is  highly 
non-linear  and  often  difficult  to  model  over  the  entire  domain  from  the  fatigue  crack 
growth  threshold  at  very  low  growth  rates  (10'9  mm/cycle)  to  the  approach  to  crack 
instability  as  the  material  toughness  is  approached  at  high  growth  rates  (near  10~2 
mm/cycle).  The  response  can  become  even  more  complex  when  aggressive  environments 
are  present  (e.g.  seawater  or  elevated  temperature)  since  under  these  conditions  the  cyclic 
frequency  and  waveform  can  influence  the  rates.  Because  of  the  general  sigmoidal  shape 
of  da/dN(AK)  several  investigators  have  employed  mathematical  function  with  a  similar 
shape  to  model  fatigue  crack  growth  rate  data  using  additional  shape  factors,  rotational 
factors  and  asymptotes  at  the  crack  growth  rate  extremes.  However,  as  has  been  pointed 
out  [22,  23]  these  asymptotic  models  are  problematic  in  the  near-threshold  regime  where 
they  often  over-estimate  the  fatigue  crack  growth  threshold  and  under  estimate  the  crack 
growth  rates,  both  of  which  are  non-conservative.  Moreover,  these  lower  asymptotes 
present  practical  problems  when  attempting  to  quantify  their  uncertainty  for  use  in 
probabilistic  analysis  of  components  and  structures. 

In  this  section,  Ti-6-4  fatigue  crack  growth  rate  data  generated  in  air  at  room  temperature 
by  three  different  laboratories  under  the  Air  Force  Office  of  Scientific  Research-High 
Cycle  Fatigue  (AFOSR-HCF)  program  [2,  4,  5]  were  analyzed.  Their  uncertainty  was  also 
quantified  such  that  the  results  could  be  used  to  characterize  fatigue  crack  growth  rates  as  a 
random  material  variable  in  a  probabilistic  life  assessment.  Near  the  fatigue  crack  growth 
threshold  accuracy  was  modeled  in  terms  of  AK,  whereas  for  higher  crack  growth  rates 
accuracy  was  modeled  in  terms  of  the  crack  growth  rates  themselves.  The  two-phase 
procedure  to  fit  the  parameters  is  outlined  below.  Due  to  the  functional  form  of  the 
multicomponent  model  employed,  the  parameters  that  describe  the  behavior  near  the 
threshold  can  be  estimated  independently  from  the  parameters  in  the  Paris  regime.  Once 
the  best-fit  parameters  are  identified,  the  variations  in  the  best-fit  parameters  among 
individual  specimens  and  from  the  different  labs  are  described  using  random  variables. 

Ideally  the  random  variable  model  is  specifically  tailored  to  the  application  and  takes 
maximum  advantage  of  the  model  form  and  physics  behind  the  data;  e.g.:  based  on 
physical  considerations,  one  may  decide  to  keep  the  Paris  slope  independent  of  the  stress- 
ratio  R.  Although  many  crack  growth  data  were  available,  rather  few  complete  fatigue 
crack  growth  rate  curves,  covering  the  entire  domain  from  the  crack  growth  threshold  to 
the  fracture  toughness  instability  limit,  were  available.  The  data  set  was  somewhat  limited, 
particularly  for  the  largest  crack  growth  rates;  nevertheless,  uncertainty  quantification 
could  be  performed  near  the  threshold  and  at  intermediate  crack  growth  rates,  which  are 
often  of  most  interest  in  fatigue  crack  growth  analyses  of  components  and  structures. 
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The  uncertainty  modeling  approach  and  its  implementation  is  described  below.  First,  the 
model  relating  fatigue  crack  growth  rate  (da/dN)  as  a  function  of  crack-tip  stress  intensity 
range  (AK)  is  described,  model  coefficients  and  their  R-dependence  are  estimated,  and 
finally  the  overall  uncertainty  is  quantified. 

4.1.1  The  Multicomponent  Model 

Typical  Ti-6A1-4V  fatigue  crack  growth  rate  data  are  shown  in  Figure  41  [2,  4,  5].  Because 
of  the  modeling  difficulties  and  errors  associated  with  asymptotic  models  a 
multicomponent  crack  growth  rate  equation  is  employed.  This  multicomponent  model  is  a 
generalization  of  the  three-component  model  previously  developed  by  Saxena  and  Hudak 
[22],  Multicomponents  are  required  for  application  to  T-6A1-4V  because  of  the  complex 
behavior  exhibited,  as  shown  in  Figure  41.  These  multicomponents  are  associated  with  the 
following  crack  growth  rate  regimes: 


Region  1 :  The  near-threshold  regime  at  very  small  crack  growth  rates  (roughly  below 
3  x  10'7  mm/cycle),  where  the  da/dN-AK  slope  is  very  steep. 

Region  2:  The  classical  Paris  crack  growth  rate  region. 

Region  3:  An  accelerated  crack  growth  region. 

Region  4:  A  crack  growth  region  similar  to  Paris  region. 

Region  5:  The  accelerated  crack  growth  corresponding  to  the  onset  of  instability. 


The  transition  between  Region  1  and  2  is  marked  by  different  physical  crack  growth 
phenomena  on  the  micro- structural  scale.  Also  the  acceleration  in  crack  growth  rates  in 
Region  3  may  be  due  to  an  environmental  effect.  Therefore,  a  series  of  two-component 
models  is  proposed:  one  for  the  Region  1-2  transition  and  one  for  the  Region  3-4  behavior. 


1  1 _ 1 _ 

C1AKm'  C2AKmi  C2  [(1  -  R)KC  ]2  AKmi~2 

r_, _ i  i  r ' 

C3AKmi  C4AKm<  C4[(l-R)Kc]2AKm 4-2 

The  model  is  set  up  such  that  the  largest  crack  growth  rate  predicted  by  the  series  of  two- 
component  models  is  selected.  Note  that  both  two-component  models  have  the  same 
fracture  roughness  Kc;  this  ensures  that  the  Region  1-2  component  produces  da/dN  rates 
below  the  Region  3-4-5  da/dN  rates.  It  should  also  be  noted  that  unless  R  is  very  high,  the 
instability  term  of  the  Region  1-2  model  is  never  exercised. 


da 

~dN 


=  max/ 
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In  this  equation,  C;  and  mi  are  constants  in  the  first  component  that  describes  crack  growth 
in  the  near-threshold  regime,  C2  and  m2  are  constants  in  the  second  component  that 
describes  crack  growth  in  the  classical  Paris  regime,  and  the  third  term  allows  for 
accelerated  growth  near  instability  due  to  the  onset  of  monotonic  failure  modes.  R  is  the 
stress  ratio,  Kmjn/Kmax,  and  Kc  is  the  fracture  toughness.  Note  that  this  three-component 
model  does  not  require  an  asymptotic  in  the  near- threshold  (Region  1),  although  as 
discussed  later  it  is  useful  to  define  an  operational  threshold  when  comparing  results  with 
data  or  defining  the  influence  of  load  ratio  in  this  regime.  Analogous  parameters  exist  for 
the  component  that  describes  the  Region  3-4-5  behavior. 


1  2  3  4  5  6  7  8  910  20  30  40  50  6070800000 
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Figure  41.  Ti-6A1-4V  Data  Obtained  at  Various  R-Values  by  Lab  1 
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4.1.2  Estimating  Model  Coefficients 

Least  squares  fits  of  da/dN(AK)  data  are  usually  accomplished  by  minimizing  the  squared 
difference  between  the  observed  and  predicted  da/dN  values.  This  is  still  applicable  in 
Region  2-3-4  and  perhaps  5;  for  Region  1,  however,  it  makes  more  sense  to  minimize  the 
error  on  AK  for  a  given  da/dN  crack  growth  rate. 

Also,  it  is  more  convenient  to  estimate  log(Q)  as  a  parameter  than  the  intercept  value  C, 
itself.  Because  the  intercept  values  are  orders  of  magnitude  different  from  the  slopes  m„ 
the  least  squares  problem  is  numerically  much  better  behaved. 

In  principle  it  is  possible,  though  numerically  somewhat  unstable,  to  estimate  all 
coefficients  simultaneously  in  a  single  least  squares  optimization  problem.  The  process 
can  be  much  improved,  however,  by  separating  the  fitting  into  the  two  physical 
components  or  regions  of  most  interest  (Region  1,  Region  2,  and  Regions  3-4-5).  Since 
the  model  is  asymptotic  near  the  upper  end  instability,  an  inequality  constraint  must  be 
introduced  for  Kc.  It  suffices  to  require  Kc  to  be  greater  or  equal  to  the  largest  measured 
value  of  AK/(1-R). 

In  Region  1,  the  error  is  measured  on  AK.  For  AK  values,  close  enough  to  the  threshold  the 
second  and  third  component  of  the  crack  growth  equation  do  not  significantly  contribute  to 
the  total  crack  growth  rate.  This  is  shown  in  Reference  [23],  Although  the  use  of  this 
approximation  (ignoring  the  second  and  third  terms  in  the  crack  growth  equation)  requires 
the  analyst  to  judge  the  crack  growth  rate  da/dN  at  which  the  transition  to  Region  2  begins, 
it  greatly  simplifies  the  computational  process. 


da/dN 


Figure  42.  Fitting  of  AK  in  Terms  of  da/dN  near  the  Threshold 
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In  Region  2,  the  error  in  measured  on  the  crack  growth  rate.  For  AK  values  in  the  Paris 
regime,  the  first  and  third  components  do  not  significantly  contribute  to  the  total  crack 
growth  rate.  Although  the  use  of  this  approximation  (ignoring  the  first  and  third  terms  in 
the  crack  growth  equation)  requires  the  analyst  the  judge  the  crack  growth  rate  da/dN  at 
which  the  transition  to  and  from  Region  2  begins,  it  greatly  simplifies  the  computational 
process.  The  quality  of  the  combined  fit  is  readily  assessed  by  plotting  the  combined  three 
components  against  the  experimental  data  (see  Figure  43).  It  should  be  noted  that  -  for 
fitting  purposes  -  there  might  be  some  gap  between  the  end  of  the  threshold  fit  region  and 
the  beginning  of  the  Paris  region.  This  “gap  region”  in  the  least-square  analysis  represents 
the  gradual  transition  between  the  threshold  and  Paris  regions. 


Figure  43.  Comparison  of  Experimental  Data  and  Model  Predictions  (Lab  1,  R  =  0.1) 
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Although  considerable  scatter  exists  between  lab  results  (crack  growth  rates  differ  by  up  to 
half  an  order  of  magnitude),  the  raw  data  for  all  labs  generally  overlap.  This  is  shown  in 
Figure  44.  Several  Lab  2  specimens  were  analyzed  (see  Figure  45).  Testing  a  specimen 
consists  of  several  start-stop  segments;  each  of  these  start-stop  sequences  introduces 
measurement  error  (see  Figure  45)  and  also  some  transient  effects  (see  Figure  50).  In  some 
cases,  it  was  readily  apparent  which  segments  should  be  eliminated  from  the  database 
(Figure  49).  In  many  cases  it  was  not  possible  to  discern  which  individual  pieces  should 
be  used  and  in  these  circumstances  all  data  were  used  and  a  single  best- fit  curve  was 
determined  (see  specimen  2  in  Figure  45).  It  should  be  noted  that  this  procedure  may 
underestimate  the  inherent  variability  of  the  data. 

Figure  43,  Figure  45,  and  Figure  46  reveal  that  the  multicomponent  model  is  capable  of 
accurately  following  the  transitions  between  all  zones  and  can  adequately  capture  the 
environmental  effects  observed  in  the  raw  data  for  R  =  0.1.  Figure  47  illustrates  the  data 
fit  at  R  =  0.5  Figure  48,  Figure  49,  and  Figure  50  illustrate  the  data  fits  for  R  =  0.8;  for  the 
high  R  very  few  data  sets  were  available  over  the  entire  da/dN(AK)  domain;  consequently, 
in  many  cases  the  fits  had  to  be  limited  to  the  Regions  1  and  2  . 

The  “bulge”  associated  with  Regions  3  and  4  is  readily  visible  in  Figure  43,  Figure  45,  and 
Figure  46  for  R  =  0. 1,  Figure  47  for  R  =  0.5  and  Figure  49  for  R  =  0.8.  The  “bulge” 
becomes  less  pronounced  with  increasing  R,  which  is  a  further  indication  that  this  is  an 
environmental  effect.  The  effect  seems  to  be  AK  related  and  not  rate  dependent;  it  occurs 
near  11-13  MPa.m0  5  for  all  R- values. 


Figure  44.  Comparison  Between  Results  from  Different  Labs 
at  R  =  0.1  (Left)  and  R  =  0.8  (Right) 
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Figure  45.  Comparison  of  Individual  Experimental  Data  Sets 
and  Model  Predictions  (Lab  2,  R  =  0.1) 
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Figure  46.  Comparison  of  Experimental  Data  and  Model  Predictions  (Lab  3,  R  =  0.1) 


53 


Figure  47.  Comparison  of  Experimental  Data  and  Model  Predictions  (Lab  1,  R  =  0.5) 
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Figure  48.  Comparison  of  Experimental  Data  and  Model  Predictions  (Lab  1,  R  =  0.8) 
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Figure  49.  Comparison  of  Experimental  Data  and  Model  Predictions  (Lab  2,  R  =  0.8) 
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Figure  50.  Comparison  of  Experimental  Data  and  Model  Predictions  (Lab  3,  R  =  0.8) 
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4.1.3  R-Dependence  of  the  Model 

Often  an  effective  AK  (AKeff)  is  used  to  correct  the  experimental  data  for  crack  closure,  in 
attempting  to  collapse  fatigue  crack  growth  rates  obtained  at  various  R- values  onto  a  single 
curve  for  which  the  crack  growth  parameters  can  then  be  estimated.  Figure  41  indicates 
that  the  necessary  closure  correction  would  not  only  be  R-dependent  but  also  very  non¬ 
linear  with  AK,  or  crack  growth  regime. 

An  alternative  approach,  which  is  proposed  here,  is  to  account  for  the  R-dependence  by 
making  certain  of  the  crack  growth  rate  modeling  parameters’  functional  dependence  on  R. 
Figure  51  shows  the  fluctuation  of  the  Paris  slope  and  log(intercept)  as  functions  of  R. 

Due  to  the  limited  data,  it  is  difficult  to  unequivocally  determine  whether  or  not  the  Paris 
slope  is  a  function  of  R.  The  following  questions  will  be  addressed: 

Are  the  mean  values  of  the  crack  growth  parameters  independent  of  R? 

Is  the  uncertainty  on  the  crack  growth  parameters  independent  of  R? 

Is  there  a  statistical  dependence  (e.g.  correlation)  between  the  crack  growth 
parameters? 

Different  modeling  approaches  are  certainly  possible.  If,  for  instance,  it  is  believed  that 
the  Paris  slope  should  be  R-independent,  the  individual  data  sets  could  be  refitted  while 
setting  the  Paris  slope  equal  to  the  average  values  obtained  for  all  specimens,  and 
subsequently  recomputing  the  intercepts  using  least  squares.  Alternatively,  based  on  the 
data  presented  in  Figure  41  one  could  require  that  the  slope  in  Region  2  and  4  be  equal 
since  both  segments  represent  the  Paris  regime.  In  that  case,  Region  3  is  interpreted  as  a 
temporary  acceleration  of  the  crack  growth  due  to  environmental  effects.  However,  none 
of  these  alternate  modeling  approaches  have  been  pursued. 
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Figure  51.  R-Dependence  of  Crack  Growth  Parameters  in  Paris  Regime 
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Figure  52  shows  the  dependence  of  the  operational  threshold  AKth  on  R.  The  operational 
threshold  is  defined  here  as  the  AK  value  at  a  crack  growth  rate  equal  to  1  x  10'7  mm/cycle 
to  be  consistent  with  the  definition  employed  by  ASTM  when  measuring  fatigue  crack 
growth  rates  and  AKth  values.  Again,  insufficient  data  are  available  to  determine  whether  a 
linear  R-dependence  can  be  applied  to  AKth. 


Figure  52.  R-Dependence  of  Threshold  Parameter  AK,h 
(Corresponding  to  da/dN  =  1  x  10‘7  mm/cycle) 


4.1.4  Uncertainty  Quantification 

It  should  be  mentioned  that  the  rapid  fluctuation  along  a  specific  da/dN  curve  is  generally 
not  of  much  interest  in  the  prediction  of  lives.  The  main  consideration  is  to  quantify  the 
uncertainty  on  the  best  fit  curves.  Although  the  data  set  is  somewhat  limited,  particularly 
when  it  comes  to  estimating  the  uncertainties  in  Regions  3-4-5,  some  uncertainty 
quantification  can  be  performed  in  Regions  1  and  2  (which  were  the  focus  of  the  HCF 
program  and  are  of  most  interest  to  fatigue  in  aircraft  components  and  structures.  A  more 
comprehensive  data  set  would  also  enable  the  uncertainties  in  Regions  3-4-5  to  be  assessed 
in  similar  fashion. 

For  Region  1  the  discussion  will  be  limited  to  the  threshold  parameter:  AKth  (Figure  52). 
The  coefficient  of  variation  on  the  threshold  parameter  AK*  coefficient  of  variation  (COV) 
is  7%  for  R  =  0. 1  and  3%  for  R  =  0.8.  Insufficient  data  are  available  to  determine  the  COV 
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of  AKth  at  R  =  0.5.  Because  of  the  additional  amount  of  closure,  and  the  micro-scale 
deformations  due  to  this  closure,  it  is  quite  typical  to  observe  more  scatter  for  R  =  0.1  than 
for  R  =  0.8.  Although  limited  data  are  available  at  both  R  =  0.1  and  R  =  0.8,  the 
standardized  normal  probability  plot  in  Figure  53  suggests  that  a  Normal  distribution  may 
be  a  reasonable  assumption  for  distribution  of  AKth.  A  summary  of  all  parameters  and 
their  uncertainty  is  given  in  Table  2. 

Table  2.  Summary  Statistics  of  the  Crack  Growth  Parameters 
in  the  Threshold  (Region  1)  and  Paris  Regimes  (Region  2) 


R  = 

0.1 

R  =  0.5 

R  = 

0.8 

Mean 

Sigma 

Value 

Mean 

Sigma 

Threshold  AKth 

4.43 

0.30 

2.87 

2.63 

0.08 

-Log-io(intercept) 

12.48 

0.57 

11.07 

11.19 

0.12 

Paris  Slope 

4.60 

0.54 

3.23 

3.87 

0.21 

In  Region  2,  the  Paris  slope  and  intercept  are  highly  correlated  as  shown  in  Figure  54. 
Specifically,  Figure  55  indicates  that  lowering  the  slope  is  offset  by  an  increase  in  the 
intercept  (i.e.  lowering  -log(intercept)).  The  correlation  coefficient  between 
-logio(intercept)  and  the  Paris  slope  is  equal  to  0.98  at  R  =  0.1,  and  equal  to  0.53  at 
R  =  0.8.  Thus,  there  is  an  almost  perfect  linear  dependence  (correlation  coefficient  =  0.98) 
between  the  Paris  slope  and  intercept  at  R  =  0.1.  This  inverse  correlation  between  the 
power-law  slope  and  intercept  has  been  observed  in  a  variety  of  materials  (for  example,  see 
Reference  [24]).  The  COV  of  the  log(intercept)  is  5%  at  R  =  0.1,  and  only  1%  at  R  =  0.8. 
The  standardized  normal  probability  plot  in  Figure  56  suggests  that  a  Normal  distribution 
may  be  a  reasonable  assumption  for  the  log(intercept)  and,  by  extension,  the  Paris  slope  as 
well  -  due  to  its  strong  correlation  with  -log(intercept). 
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Figure  53.  Standardized  Normal  Probability  Plot  for  the  Threshold  AKth 
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Figure  54.  Paris  Coefficients  for  Various  R-Values 
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The  plot  indicates  very  strong  correlation  between  log(intercept)  and  Paris  Slope  at  R  =  0.1 . 


Figure  55.  Joint  Plot  of  Standardized  Intercept  and  Slope  in  Region  2  (Paris  Region) 


Figure  56.  Standardized  Normal  Probability  Plot  for  Coefficients  in  Region  2 

(Paris  Region) 
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4.2  Statistical  Modeling  of  Mission  Variability 

The  crack  growth  life  associated  with  a  gas  turbine  engine  disk  is  influenced  by  a  number 
of  random  variables  (e.g.,  material  properties,  inherent/induced  anomalies,  applied  stress, 
among  others).  The  dispersion  associated  with  each  of  these  variables  influences  the 
probability  of  fracture.  The  accuracy  of  predicted  risk  is  highly  dependent  on  the  statistical 
models  used  to  represent  each  of  the  individual  random  variables.  Parametric  models  (i.e., 
probability  density  functions  such  as  normal,  lognormal,  Weibull)  are  commonly  used  to 
represent  the  random  variables  in  which  the  accuracy  of  each  variable  is  dependent  on  a 
regression  fit  of  the  data.  Since  experimental  data  are  typically  most  plentiful  near  the 
mean  value,  the  regression  curve  can  be  biased  towards  the  mean  and  is  often  a  poor  fit 
near  the  tails  of  the  distribution.  Also,  parametric  models  contain  values  that  are  often 
physically  impossible  to  achieve  in  physical  systems  (e.g.,  ±  co).  In  addition,  although  the 
regression  approach  can  be  used  to  estimate  the  degree  of  correlation  among  two  or  more 
random  variables,  it  may  not  adequately  address  general  relationships  among  multiple 
random  variables  (e.g.,  multimodality,  skewed  discontinuities,  among  others). 

In  this  section,  a  nonparametric  method  is  presented  for  density  estimation  of  rainflow 
stress  profiles  associated  with  aircraft  gas  turbine  engine  usages.  A  brief  summary  of  the 
adaptive  kernel  density  estimator  is  presented  including  details  regarding  a  specific 
algorithm.  Algorithm  performance  is  illustrated  using  data  generated  from  parametric  uni- 
and  multi-variate  probability  density  functions  and  is  compared  to  parametric  estimates 
based  on  the  same  data.  The  adaptive  kernel  density  estimator  is  applied  to  rainflow  stress 
pairs  associated  with  F-16/F100  engine  usages  and  compared  to  parametric  estimates, 
including  a  hybrid  approach  based  on  separate  treatment  of  maximum  stress  pairs 
associated  with  engine  startup  and  aircraft  takeoff. 

Usage  histories  for  commercial  aircraft  are  often  well  established  with  known  stress  pair 
magnitudes.  A  stress  scatter  factor  can  be  used  to  account  for  the  bias  between  predicted 
and  measured  stress  values.  This  factor  can  be  modeled  as  a  random  variable  to  address 
the  variability  associated  with  the  measured  values.  As  shown  in  Figure  57(a),  the 
influence  of  stress  scatter  on  the  probability  of  fracture  is  commonly  estimated  by  applying 
it  repeatedly  to  the  deterministic  stress  values  over  the  entire  life  of  a  disk. 

For  military  usages,  the  number  and  magnitude  of  stress  pair  values  are  not  always  well 
defined,  and  may  have  greater  variability  compared  to  commercial  usages.  Since  the  stress 
pair  values  could  vary  significantly  from  flight  to  flight,  a  more  accurate  model 
(Figure57(b))  is  needed  to  estimate  the  influence  of  stress  variability  on  the  probability  of 
fracture. 

For  gas  turbine  engines,  the  applied  stress  is  a  stochastic  variable  that  is  often  difficult  to 
measure  and  predict.  Two  primary  approaches  have  been  proposed  for  density  estimation 
of  rainflow  stress  values:  parametric  and  nonparametric.  The  parametric  approach  consists 
of  a  regression  fit  of  the  data  plotted  on  a  cumulative  probability  scale  [25,  26,  27].  The 
primary  drawback  to  this  approach  is  that  it  is  often  difficult  to  apply  to  multimodal 
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densities.  On  the  other  hand,  nonparametric  methods  provide  a  direct  fit  to  the  data  and  are 
routinely  used  to  estimate  multimodal  densities.  The  kernel  density  estimator,  the  most 
commonly  used  nonparametric  method,  has  been  recently  applied  to  rainflow  stress  data 


Figure  57.  Approaches  for  Modeling  Variability  Associated  with  Applied  Loads: 
(a)  Stress  Bias  Model,  and  (b)  Cycle-by-Cycle  Usage  Variability  Model 


A  number  of  nonparametric  methods  can  be  used  to  estimate  the  probability  density 
associated  with  a  data  set,  such  as  the  nearest  neighbor  method  [30],  convex  intensity 
[31],  and  orthogonal  series  estimators  [32],  among  many  others.  Many  of  the 
nonparametric  methods  are  related  to  the  histogram,  which  is  commonly  used  to  provide  a 
measure  of  the  relative  occurrences  of  discrete  data  values.  However,  as  shown  in 
Figure  58,  histograms  are  highly  dependent  on  the  bin  width  selected  for  sampling. 
Furthermore,  a  histogram  is  not  a  true  probability  density  function,  and  therefore  cannot  be 
used  for  probabilistic  computations. 
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Figure  58.  Bin-Width  Dependency  Associated  with  Histogram  Density  Estimates: 
(a)  Bin  Width  =  5  Units,  and  (b)  Bin  Width  =  10  Units 


The  kernel  density  estimator  is  one  of  the  most  commonly  used  nonparametric  methods  to 
estimate  probability  densities.  As  illustrated  in  Figure  59(a),  each  data  point  is  assumed  to 
have  an  associated  probability  density  within  a  specified  range  of  values.  If  data  points  are 
relatively  close  together,  their  probability  densities  can  overlap  (Figure  59(b)).  The 
probability  density  at  a  given  location  can  be  estimated  as  the  sum  of  the  contributions  of 
the  kernels  associated  with  each  data  point  [30]: 
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where  tj,  t2  =  coordinates  of  evaluation  point,  n  =  number  of  data  points,  h  =  window 
width  (smoothing  parameter),  K  =  kernel  estimator,  and  Xu,  Xtj  =  coordinates  of  data 
point  i. 
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Equation  13  provides  an  initial  density  estimate  associated  with  a  fixed  window  width  h. 
The  window  width  can  be  adaptively  adjusted  to  account  for  the  number  of  data  points 
within  a  region  using  the  following  equation  [32]: 
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Figure  59.  Nonparametric  Approach  to  Probability  Density  Estimation: 
(a)  Epanechnikov  Kernal  Density  at  Individual  Data  Points,  and 
(b)  Combined  Density  of  All  Data  Points 
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where  l  is  a  bandwidth  parameter  that  is  identified  adaptively  based  on  previous 
probabilistic  density  function  (PDF)  estimates: 
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An  Epanechnikov  kernel  function  can  be  used  to  describe  the  probability  density 
associated  with  rainflow  stress  values  [29,  32]: 
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An  adaptive  kernel  numerical  algorithm  based  on  Equations  13-16  was  developed  and 
verified  using  a  number  of  established  probability  density  functions.  For  example, 
consider  the  univariate  lognormal  distribution  (LN(3,2))  shown  in  Figure  60(a).  The  initial 
kernel  PDF  estimate  based  on  Equation  13  (and  200  Monte  Carlo  samples  from  the  parent 
distribution)  is  very  similar  to  the  parent  PDF  near  the  median,  but  differs  near  the  right 
tail  of  the  distribution.  However,  when  the  window  width  is  adaptively  adjusted  using 
Equation  14,  it  converges  to  the  parent  distribution  over  most  of  the  right  tail.  The 
parametric  density  estimate  (LN(2.89,  1.81)  is  very  similar  to  the  parent  both  at  the  mode 
and  tails  of  the  distribution. 

The  algorithm  was  applied  to  the  bivariate  lognormal  distribution  (LN(3,1.5)  +  LN(9,1.5)) 
shown  in  Figure  60(b).  Both  the  initial  and  adaptive  kernel  density  estimates  predict  the 
location  of  the  modes  and  provide  reasonable  treatment  at  the  tails  of  the  distribution.  On 
the  other  hand,  the  parametric  estimate  is  a  poor  fit  both  at  the  modes  and  tails  of  the 
distribution. 

The  algorithm  was  also  verified  for  the  bivariate  normal  distribution  shown  in  Figure  61 
(X  =  N(100,10),  7=  N(  100, 10),  uncorrelated  X  and  T)  based  on  200  samples  from  the 
parent  distribution.  The  parametric  fit,  shown  in  Figure  61(a),  is  nearly  identical  to  the 
parent  distribution.  The  kernel  density  estimate  (Figure  61(b))  has  an  overall  shape  that  is 
similar  to  the  parent,  but  is  quite  different  in  the  tails  of  the  distribution.  Standard 
algorithms  are  available  to  assess  the  accuracy  of  the  kernel  estimate  [30]  which  should  be 
strongly  considered  before  applying  the  density  estimate  for  risk  predictions. 

The  adaptive  kernel  algorithm  was  applied  to  the  density  estimation  of  rainflow  stress  pairs 
associated  with  actual  aircraft  gas  turbine  engine  usages.  Raw  RPM  values  were  obtained 
from  engine  flight  data  recorder  data  [33]  for  a  variety  of  F-16/F100  missions.  RPM 
values  are  shown  in  Figure  62  for  representative  flights  of  missions  classified  as  “Live 
Fire”  and  “Instruments  &  Navigation”.  The  RPM  values  were  converted  to  stress  values 
using  an  empirical  algorithm  including  a  shakedown  procedure  to  address  stress  values  that 
fall  above  the  material  yield  stress  [34],  The  stress  values  were  sorted  into  min-max  pairs 
using  an  established  rainflow  cycle  counting  algorithm  [35]. 

As  shown  in  Figure  62,  a  typical  engine  usage  history  can  consist  of  thousands  of  RPM 
data  points.  However,  the  data  can  often  be  reduced  to  a  few  hundred  values  following 
application  of  the  rainflow  cycle  counting  algorithm.  For  example,  the  rainflow  stress 
pairs  associated  with  each  of  the  individual  flights  shown  in  Figure  63  are  based  on  usages 
similar  to  the  ones  shown  in  Figure  62.  Most  of  the  stress  pairs  are  clustered  in  a  region  of 
minimum  and  maximum  stress  values  ranging  from  350  to  750  MPa.  However,  a  few 
stress  pairs  associated  with  the  aircraft  major  stress  cycle  are  located  in  the  region  where 
the  minimum  stress  is  zero  (associated  with  engine  startup  and  aircraft  takeoff). 
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(a) 


(b) 


Figure  60.  Comparison  of  Parametric  and  Nonparametric  Density  Estimation 
Methods  for  (a)  Univariate  Unimodal  Lognormal  Parent  Distribution, 
and  (b)  Univariate  Bimodal  Parent  Distribution 
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Figure  61.  Comparison  of  Density  Estimation  Methods  for  a  Multivariate  Normal 
Distribution:  (a)  Parametric  Density  Estimate,  and  (b)  Kernel  Density  Estimate 
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Figure  62.  F-16/F100  Aircraft  Engine  Usage  Histories  Associated  with  Representative 
Flights  of  the  Live  Fire  and  Instruments  &  Navigation  Missions 
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Figure  63.  Rainflow  Stress  Pairs  Associated  with  Representative  Flights  for  the 
F-16/F100  Aircraft  Live  Fire  Usage  Classification 


If  the  stress  pairs  can  be  completely  described  by  a  probability  density  function,  then  the 
influence  of  usage  variability  on  the  probability  of  fracture  can  be  quantified  by  sampling 
directly  from  the  rainflow  stress  PDF.  The  desired  PDF  must  provide  a  probability  density 
for  each  region  of  the  population  that  is  consistent  with  the  clustering  of  data  points  shown 
in  Figure  63.  It  must  attempt  to  describe  the  population  from  which  the  data  points  were 
drawn  to  reflect  the  influence  of  data  points  not  included  in  the  current  sample. 

The  stress  pairs  associated  with  startup/takeoff  can  be  considered  separately  from  the 
remaining  data  because  they  must  be  applied  only  once  per  flight  (this  cannot  be 
guaranteed  if  they  are  included  in  the  sampling  population).  In  addition,  since  the 
minimum  stress  values  associated  with  these  stress  pairs  are  always  zero  (to  reflect  the  zero 
RPM  value  prior  to  engine  startup),  they  can  be  modeled  as  a  univariate  (and  probably 
unimodal)  distribution.  A  parametric  probability  density  estimate  of  the  sub-population  of 
Live  Fire  rainflow  stress  pairs  that  does  not  include  values  associated  with  startup/takeoff 
is  shown  in  Figures  64(a)  and  (c).  The  parametric  PDF  appears  to  address  the  correlation 
structure  of  minimum  and  maximum  stress  values  presented  in  Figure  63. 
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A  kernel  density  estimate  of  the  Live  Fire  mission  rainflow  stress  pairs  is  shown  in 
Figures  64(b)  and  (d).  It  provides  explicit  treatment  of  the  stress  pairs  associated  with 
engine  startup,  as  indicated  by  the  stress  values  that  are  separated  from  the  primary  cluster 
that  appear  near  zero  on  the  min  stress  axis.  It  also  appears  to  capture  the  correlation 
structure  of  the  minimum  and  maximum  stress  values.  It  is  interesting  to  note  that  the 
parametric  method  predicts  a  unimodal  density  for  the  rainflow  data,  whereas  the  kernel 
method  predicts  a  bimodal  density.  Since  it  is  expected  that  extreme  stress  values  will 
have  more  influence  on  the  probability  of  fracture,  it  may  be  more  appropriate  to  assess  the 
accuracy  of  the  methods  in  the  tails  of  the  distributions. 

In  Figure  63,  it  can  be  observed  that  the  maximum  and  minimum  stress  pair  values 
associated  with  the  Live  Fire  mission  appear  to  be  correlated  random  variables.  If  these 
values  were  modeled  as  random  variables  for  risk  assessment,  a  correlated  sampling 
technique  would  be  required  to  capture  this  relationship.  Although  correlated  sampling  is 
well  established  for  parametric  random  variables,  significant  errors  may  be  introduced 
when  the  random  variables  are  non-Gaussian.  Correlated  sampling  is  not  as  mature  for 
nonparametric  variables,  and  the  errors  and  limitations  are  not  as  well  known. 

Common  correlated  sampling  techniques  involve  transformation  of  a  vector  of  correlated 
variables  into  a  vector  of  uncorrelated  variables  using  matrix  algebra  methods  (e.g., 
Rosenblatt  transformation)  that  are  similar  to  coordinate  transformation  [36],  In  some 
instances,  a  coordinate  transformation  can  be  identified  for  the  original  random  variables 
that  eliminates  the  need  for  uncorrelated  sampling.  For  example,  consider  the  max-min 
stress  pair  values  shown  in  Figure  63.  If  the  coordinate  system  is  transformed  from  max- 
min  stress  to  delta  stress-R  ratio  (Figure  65),  it  can  be  observed  that  the  data  exhibit  strong 
negative  correlation  (correlation  coefficient  =  -0.992).  If  it  is  assumed  that  the  transformed 
variables  are  perfectly  negatively  correlated,  then  uncorrelated  sampling  techniques  can  be 
used  to  simulate  these  variables.  This  approach  avoids  the  errors  associated  with 
correlated  sampling,  and  applies  to  both  parametric  and  nonparametric  variables.  The 
density  estimation  of  stress  pair  values  is  greatly  simplified  when  modeled  as  delta  stress 
and  R  ratio.  If  the  PDF  is  simulated  using  uncorrelated  sampling  techniques,  the  density 
estimation  can  be  performed  on  each  random  variable  individually. 
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Figure  64.  Comparison  of  Density  Estimation  Methods  for  Rainflow  Stress  Pairs 
Associated  with  Live  Fire  Usage  Histories:  (a)  and  (c)  Parametric  Density  Estimate 
(No  Engine  Startup  Stress  Pairs);  (b)  and  (d)  Kernel  Density  Estimate 
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Figure  65.  The  Rainflow  Stress  Values  for  the  Live  Fire  Mission  Exhibit  Strong 
Negative  Correlation  When  Transformed  to  Delta  Stress  Versus  R-Ratio  Space 


When  sampling  transformed  random  variables,  special  care  must  be  taken  to  ensure  that 
the  stress  values  selected  are  consistent  with  the  field  data  values.  For  example,  consider 
the  delta  stress  versus  stress  pair  number  values  for  the  Live  Fire  mission  flights  shown  in 
Figure  66.  The  rainflow  procedure  provides  stress  pair  values  that  decrease  in  magnitude 
with  increasing  stress  pair  number  for  an  individual  flight.  The  range  of  stress  pair  values 
associated  with  each  stress  pair  number  must  also  be  considered.  In  Figure  67,  mean,  5th 
and  95th  percentile  delta  stress  magnitudes  are  shown  for  each  stress  pair  number  for  the 
Live  Fire  mission.  Based  on  the  data  shown  in  Figure  67,  it  appears  that  the  delta  stress 
values  associated  with  each  stress  pair  could  be  modeled  as  random  variables. 

The  main  descriptors  (mean,  standard  deviation,  COV)  for  each  stress  pair  are  shown  in 
Figure  68  (shown  as  normalized  values).  The  relationship  among  delta  stress  values 
associated  with  successive  stress  pairs  must  also  be  considered.  In  Figure  69,  delta  stress 
values  are  shown  for  two  representative  stress  pairs  of  the  Live  Fire  mission,  indicating 
strong  correlation  among  delta  stress  values.  This  relationship  is  confirmed  in  Figure  70, 
where  delta  stress  values  are  shown  for  representative  individual  flights.  It  can  be 
observed  that  all  of  the  delta  stress  values  are  either  above  the  mean  or  below  the  mean 
stress  values  for  at  least  some  of  the  flights.  However,  it  should  be  noted  that  for  some  of 
the  flights,  the  delta  stress  values  appear  both  above  and  below  the  mean  values,  which 
suggests  that  stress  pair  values  are  partially  correlated. 
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Figure  66.  Delta  Stress  and  Stress  Pair  Number  Associated  with 
Individual  Flights  of  the  Live  Fire  Mission 


Stress  Pair 


Figure  67.  Delta  Stress  Versus  Stress  Pair  Number  for  All  Flights 
Associated  with  the  Live  Fire  Mission 
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Figure  68.  Normalized  Mean,  Standard  Deviation,  and  Coefficient  of  Variation 
(COV)  Associated  with  Delta  Stress  Values  for  the  Live  Fire  Mission 
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Figure  69.  Delta  Stress  Values  Associated  with  Two  Stress  Pair  Numbers 
for  Each  Flight  of  the  Live  Fire  Mission 
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Note  that  all  of  the  delta  stress  values  for  an  individual  flight  may 
appear  above  or  below  the  mean  mission  values. 


Figure  70.  Delta  Stress  Values  for  Selected  Flights  of  the  Live  Fire  Mission 
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4.3  Probabilistic  Mission  Identification  (PMI) 

Fatigue  life  prediction  of  aircraft  gas  turbine  engine  rotating  components  requires  estimates 
of  the  applied  stress  values  throughout  the  life  of  the  component.  These  values  may  vary 
considerably  from  flight-to-flight,  and  are  highly  dependent  upon  the  mission  type. 
Unfortunately,  engine  flight  data  recorders  currently  do  not  have  the  capability  to  identify 
the  mission  type  for  a  given  flight  history.  Instead,  the  aircraft  pilot  must  manually  enter 
the  mission  type  in  a  log  book,  which  can  lead  to  errors  and  omissions  of  this  critical 
information.  In  addition,  a  single  flight  may  have  characteristics  that  are  similar  to 
multiple  mission  types.  Thus,  the  mission  identification  problem  consists  of  selecting  the 
most  likely  mission  type  from  a  set  of  established  standard  missions.  An  automated 
mission  identification  method  would  greatly  improve  risk  assessment  estimates, 
particularly  if  it  were  based  on  quantitative  methods. 

A  quantitative  probabilistic  method  has  been  developed  to  identify  the  mission  type  for  a 
given  flight  history.  The  method  is  based  on  a  comparison  of  the  rainflow  stress  PDF  of  a 
specific  flight  history  with  previously  developed  PDFs  of  standard  mission  types.  The 
method,  called  PMI,  predicts  the  most  likely  mission  type  based  on  volume  integration  of 
the  region  common  to  both  the  flight  history  and  the  standard  mission. 

A  description  of  the  probability  densities  associated  with  one  or  more  random  variables  is 
necessary  for  probabilistic  risk  assessment.  A  number  of  parametric  models  are  available 
to  estimate  probability  density  based  on  a  statistical  fit  of  the  available  data.  The 
appropriate  distribution  is  based  on  the  statistical  error  associated  with  the  data  fit  that  is 
quantified  using  classical  goodness-of-fit  tests  (e.g.,  chi-square,  Kolmogorov- Smirnov, 
among  others  [25,  26]).  These  tests  provide  a  quantitative  measure  of  the  relationship 
among  a  data  set  and  a  parametric  distribution. 

The  mission  identification  problem  consists  of  identifying  the  distribution  that  most  closely 
resembles  the  available  data.  Consider  mission  Mt  (ith  unknown  mission)  and  fleet  Fj  (/th 
fleet  of  known  mission  type)  shown  in  Figure  71.  The  volume  of  the  region  over  which 
the  probability  densities  of  the  stress  values  (s„  S2)  associated  with  mission  Mi  and  fleet  Fj 
overlap  can  be  used  as  a  measure  of  the  likelihood  that  mission  Mi  is  a  member  of  fleet  Fj. 
The  overlap  probability  can  be  expressed  as: 
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Noting  that  the  denominator  in  Equation  17  has  a  value  of  1  when  integrated  over  the 
entire  domain  of  Fj,  Equation  17  can  be  simplified  as: 
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Note  that  Equation  18  applies  to  both  parametric  and  nonparametric  probability  densities. 


Figure  71.  The  Volume  of  the  Region  Over  Which  the  Probability  Densities  of  the 
Stress  Values  (Sj,  .v2)  Associated  with  Mission  M,-  and  Fleet  Fj  Overlap  Can  Be  Used 
as  a  Measure  of  the  Likelihood  That  Mission  M,-  is  a  Member  of  Fleet  Fj 


The  PMI  method  was  used  to  identify  the  mission  type  associated  with  actual  military 
aircraft  engine  usages.  Raw  RPM  values  (Figure  63)  were  obtained  from  engine  flight 
recorder  data  and  converted  to  stress  values  using  an  algorithm  that  accounts  for  local 
yielding  [37,  38].  An  established  rainflow  counting  technique  [35]  was  used  to  sort  the 
stress  data  into  min-max  pairs.  Rainflow  stress  values  are  shown  in  Figure  72  for  two 
representative  mission  types  (Live  Fire,  Instruments  and  Navigation).  The  grouping  of 
data  is  similar  for  both  mission  types,  with  most  values  located  within  a  region  of 
minimum  and  maximum  stress  range  values  ranging  from  360  to  720  MPa. 


Probability  densities  were  estimated  using  the  adaptive  kernel  method  based  on  10  flights 
for  each  mission  type,  shown  in  Figure  73.  The  probability  density  estimates  for  both 
mission  types  appear  to  capture  the  correlation  among  minimum  and  maximum  stress 
values  associated  with  the  original  data  set.  The  PDF  associated  with  the  Live  Fire 
mission  (Figure  73(a))  appears  to  have  a  single  dominant  mode,  whereas  the  PDF 
associated  with  the  Instruments  &  Navigations  mission  (Figure  73(b))  is  characterized  by 
two  modes  with  similar  intensities. 


To  evaluate  the  performance  of  the  PMI  method,  rainflow  stress  probability  densities 
associated  with  usages  of  known  mission  type  were  compared  to  fleets  of  usages  of  known 
mission  type.  Two  test  missions  of  known  type  (Live  Fire,  Instruments  and  Navigation) 
were  selected  to  evaluate  the  PMI  method.  Probability  densities  obtained  using  the 
adaptive  kernel  method  are  shown  in  Figure  74,  and  probability  contours  associated  with 
the  test  missions  are  shown  in  Figure  75. 
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Figure  72.  Rainflow  Stress  Pair  Values  Associated  with  the  Following  Aircraft 
Usage  Classifications:  (a)  Live  Fire,  and  (b)  Instruments  &  Navigation 
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Figure  73.  Estimated  Joint  Probability  Densities  Associated  with  the  Following 
Aircraft  Usage  Classifications:  (a)  Live  Fire,  and  (b)  Instruments  &  Navigation 
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Figure  74.  Probability  Densities  Associated  with  Individual  Flights  of  Known  Mission 
Type  Were  Used  to  Illustrate  the  Performance  of  the  Probabilistic  Mission 
Identification  (PMI)  Method:  (a)  Live  Fire  Test  Mission  1,  and 
(b)  Instruments  &  Navigation  Test  Mission  1 
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Figure  75.  Probability  Contours  Associated  with  Individual  Flights  of  Known  Mission 
Type:  (a)  Live  Fire  Test  Mission  1,  and  (b)  Instruments  &  Navigation  Test  Mission  1 
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The  PMI  method  was  first  applied  to  a  test  mission  of  type  Live  Fire.  Probability  densities 
and  contours  associated  with  the  Live  Fire  test  mission  and  the  Live  Fire  and  Instruments 
&  Navigation  fleets  are  shown  in  Figures  76(a)  and  77(a),  respectively.  The  bimodal 
shape  of  the  Live  Fire  test  mission  suggests  that  it  may  belong  to  the  Instruments  and 
Navigation  fleet.  However,  when  the  probability  densities  of  the  overlap  among  the  test 
mission  and  the  fleets  are  considered  (Figures  77(a)  and  (b)),  the  region  of  relatively  high 
density  shown  in  Figure  77(a)  suggests  that  the  test  mission  belongs  to  the  Live  Fire  fleet. 
Volume  integration  of  the  overlap  density  (Table  3)  confirms  that  the  test  mission  most 
likely  belongs  to  the  Live  Fire  fleet. 

Table  3.  Overlap  Probabilities  Associated  with  Test  Missions 


Test  Mission 

Overlap  Probability 

Live  Fire 

Instr.  & 

Fleet 

Navigation 

Fleet 

Live  Fire 

0.254 

0.193 

Instruments  &  Navigation 

0.057 

0.073 

The  PMI  method  was  next  applied  to  a  test  mission  of  type  Instruments  and  Navigation.  In 
contrast  with  the  Live  Fire  test  mission,  the  Instruments  and  Navigation  test  mission  has  a 
unimodal  appearance  which  suggests  that  it  may  belong  to  the  Live  Fire  fleet.  However,  as 
shown  in  Figure  78,  the  peak  density  of  the  test  mission  is  very  close  to  one  of  the  modes 
of  the  Instruments  and  Navigation  fleet.  Volume  integration  of  the  overlap  density 
(Figure  79,  Table  3)  confirms  that  the  test  mission  most  likely  belongs  to  Instruments  and 
Navigation  fleet.  For  this  example,  the  PMI  method  correctly  identified  the  mission  type 
associated  with  both  test  missions. 

PMI  has  a  critical  role  in  the  probabilistic  treatment  of  aircraft  engine  usage  as  indicated  in 
the  conceptual  framework  shown  in  Figure  80.  The  framework  consists  of  the  following 
four  stages:  (1)  data  retrieval,  (2)  mission  identification,  (3)  stress  characterization,  and 
(4)  risk  prediction.  In  the  data  retrieval  stage,  raw  RPM  values  from  engine  sensors  are 
converted  to  stress  values  and  sorted  using  a  standard  rainflow  procedure.  In  the 
identification  stage,  PMI  is  used  to  predict  the  most  likely  mission  type  of  each  flight 
based  on  volume  integration  of  rainflowed  stress  pair  values  associated  with  the  flight 
history  and  a  standard  mission.  In  the  characterization  stage,  probability  densities 
associated  with  each  mission  type  are  quantified  using  an  adaptive  kernel  approach  that 
addresses  extreme  behavior  near  the  tails  of  the  distributions,  and  can  easily  be  updated  as 
new  data  are  acquired.  Stress  data  are  transformed  to  a  stress  ratio  versus  delta  stress 
coordinate  system  to  provide  improved  numerical  treatment  of  the  correlation  among 
rainflowed  stress  pair  values.  In  the  final  stage  (prediction),  future  usage  is  simulated 
based  on  values  obtained  from  the  mission  probability  density  functions  and  future  mission 
mix  values  from  fleet  management  projections. 
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Live  Fire  Fleet  PDF 


Figure  76.  Probability  Contours  Associated  with  an  Individual  Flight  of  Known 
Mission  Type  (Live  Fire  Test  Mission  1)  and  (a)  Live  Fire  Mission,  and 
(b)  Instruments  &  Navigation  Mission 
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Figure  77.  Overlap  Probabilities  for  Live  Fire  Test  Mission  1  and  the  Live  Fire 
and  Instruments  &  Navigation  Missions,  (a)  and  (b),  Respectively 
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Figure  78.  Probability  Contours  Associated  with  an  Individual  Flight  of  Known 
Mission  Type  (Instruments  &  Navigation  Test  Mission  1)  and  (a)  Live  Fire 
Mission,  and  (b)  Instruments  &  Navigation  Mission 
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Figure  79.  Overlap  Probabilities  for  Instruments  &  Navigation  Test  Mission  1  and 
the  Live  Fire  and  Instruments  &  Navigation  Missions,  (a)  and  (b),  Respectively 
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Figure  80.  A  Conceptual  Framework  for  Probabilistic  Treatment 
of  Aircraft  Engine  Usage 


4.4  Probabilistic  Sensitivity  Methods 

The  development  of  probabilistic  sensitivities  is  frequently  considered  an  essential 
component  of  a  probabilistic  analysis  and  often  critical  towards  understanding  the  physical 
mechanisms  underlying  failure  and  modifying  the  design  to  mitigate  and  manage  risk. 

One  useful  sensitivity  is  the  partial  derivative  of  the  probability-of- failure  with  respect  to 
the  parameters  of  the  input  random  variables’  PDF  e.g.,  dP,  I  dpi,  dPf  Ida,  where  pi  is  the 

mean  value  and  a  is  the  standard  deviation.  The  score  function  method,  described  below, 
is  a  convenient  method  to  compute  these  sensitivities  because  the  equations  are  formulated 
such  that  the  same  Monte  Carlo  samples  used  to  compute  the  probability-of-failure  of  the 
system  can  also  be  used  to  compute  the  probabilistic  sensitivities.  Therefore,  the 
sensitivities  are  obtained  for  negligible  additional  cost. 

4.4.1  Methodology 

The  score  function  method  has  been  published  previously  [39]  but  was  extended  in  this 
program  to  consider  issues  pertinent  to  DARWIN  analyses;  namely,  the  zone-based 
analysis  approach. 

The  essence  of  the  method  for  independent  random  variables  and  a  single  zone  is  presented 
below.  More  details  on  the  DARWIN  implementation  are  contained  in  the  appendix  in 
Ref.  [40], 

The  probability-of-failure  is  computed  as 

Pf  =  I  /  (x)  •  dx  (19) 

j  Jg(\)<0Jxy  '  v  ' 
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where  X  is  a  vector  of  random  variables,  /X(x)is  the  joint  density  function  of  X,  and  g 
is  the  limit  state  function  defined  such  that  g(x)  <  0  denotes  failure. 

Inserting  the  indicator  function 7(x),  which  is  defined  as  equal  to  one  if  g(x)  <  0  and  zero 
otherwise, 


Pf=  J  f(x)/x(x)Jx 


(20) 


The  derivative  of  Pf  with  respect  to  a  parameter  of  the  joint  PDF  e.g.,  mean,  standard 
deviation,  correlation  coefficient,  shape,  scale  factor,  etc.,  denoted  0,  can  be  written 

=  J  7(x)^y  ,1  /X(x)^x  =  ^[/(X)-K^x)]  (21) 

30  d0  /X(x) 


where 


g/x(x)  1  ^n(/x(x)) 

dd  /X(x)  at? 


denotes  the  kernel  functions. 

As  examples,  for  independent  random  variables,  the  kernel  functions  for  a  normal  random 
variable  are 


where 
X  and 


1 

<J2  —  (7 

Ka=-{U2- 1) 


(23) 


<7 

U=(x-ju)/<t,  and  for  a  lognormal  random  variable  defined  in  terms  of  the  median 
coefficient  of  variation  COY, 


ln(x)  -  ln(X) 
lln(l  +  C0F2) 

COV  •  ^ln(l  +  COVry  (ln(!)  -  ln(.r))2  J 
(l  +  COV2  M  +  COF2)2 


(24) 
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4.4.2  DARWIN  Implementation 

A  DARWIN  hard  alpha  analysis  consists  of  three  stages:  single  zone  conditional  risk 
assessment  (assumes  a  defect  is  present  in  the  zone),  a  single  zone  unconditional  risk 
assessment  (considers  the  probability  of  having  a  defect),  and  a  multi-zone  risk 
assessment.  The  sensitivity  results  from  each  level  feed  into  the  next  DARWIN  analysis 
for  hard  alpha  has  three  random  variables:  exceedance  curve,  stress  scatter  and  life 
propagation  scatter.  The  mathematical  development  of  the  requisite  equations  are  provided 
in  the  appendix  in  Ref.  [40]  (particularly  Equations  7,  8,  and  15  in  the  referenced 
document.) 

Stress  and  propagation  scatter  are  represented  as  lognormal  random  variables;  therefore, 
sensitivities  with  respect  to  the  median  and  COV  are  computed.  The  exceedance  curve  is 
nonparameteric  tabular  input,  defined  in  part  by  the  min  and  max  bounds  ,  amm  and  anydx , 
of  the  exceedance  curve.  Sensitivities  with  respect  to  these  bounds  are  computed  by 
DARWIN.  The  bounds  of  the  exceedance  curve  (e.g.,  the  standard  hard  alpha  anomaly 
distributions  included  in  the  Federal  Aviation  Administration  Advisory  Circular  33.14-1 
[41])  have  generally  been  set  with  the  expectation  that  they  are  far  enough  above  or  below 
the  range  of  expected  behavior  as  to  have  no  significant  impact  on  the  calculation. 
Therefore,  it  is  anticipated  that  the  sensitivities  should  be  insignificant;  otherwise,  the 
exceedance  curve  should  be  extended  as  needed  to  reduce  the  sensitivity. 

Since  the  sensitivities  are  computed  using  Monte  Carlo  sampling,  the  sensitivity  estimates 
are  random  variables.  An  estimate  of  the  variance  of  the  sensitivities  were  derived.  In 
addition,  the  distribution  of  the  sensitivities  follows  a  normal  distribution.  As  such,  the 
95%  lower  and  upper  confidence  bounds  were  developed  and  implemented.  That  is,  it  is 
95%  likely  that  the  true  sensitivity  falls  with  the  values  computed  by  DARWIN.  The 
bounds  are  set  as 


S0  - 1 .96Vf  <  Se  <  S0  + 1 ,96VC 


(25) 


where  V  is  the  variance  estimate.  Equations  for  V  are  provided  in  the  Ref.  [40]. 

4.4.3  Numerical  Example 

The  model  consists  of  a  rectangular  cross-section  titanium  ring  modeled  by  24  zones 
(Figure  81).  The  titanium  ring  is  under  rotation,  and  has  an  external  load  applied  at  the  rim 
to  simulate  blade  loading.  The  maximum  principle  stress  is  in  the  circumferential 
direction.  The  crack  grows  under  constant  amplitude  cyclical  loading  e.g.,  start  up-shut 
down,  until  fracture,  defined  as  when  the  stress  intensity  factor  exceeds  the  fracture 
toughness.  The  probability-of-fracture  of  the  disk  before  twenty  thousand  cycles  is 
computed  using  Monte  Carlo  sampling. 
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The  sensitivity  results  for  a  single  zone  unconditional  probability  of  fracture  (POF)  are 
compared  against  finite  difference  estimates  in  Table  4.  The  finite  difference  estimates  are 
computed  by  perturbing  each  parameter  by  a  small  percentage  and  rerunning.  We  note  that 
it  is  exceedingly  laborious  to  obtain  finite  difference  estimates  in  that  a  very  large  number 
of  Monte  Carlo  samples  must  be  used  in  order  to  prevent  the  difference  in  the  POF  from 
being  masked  by  sampling  variance,  or  a  very  refined  importance  sampling  method  must 
be  used.  In  addition,  a  separate  analysis  must  be  performed  for  each  sensitivity  estimate.  In 
the  examples  below,  six  additional  analyses  are  required  to  obtain  the  finite  difference 
estimates,  whereas,  the  sensitivity  equations  provide  all  sensitivities  in  a  single  analysis 
without  any  additional  sampling. 

The  sensitivity  results  for  the  disk  (with  24  zones)  are  compared  against  finite  difference 
estimates  in  Table  5.  In  both  tables,  the  comparison  between  the  sensitivity  equations  and 
finite  difference  estimates  are  quite  good. 

Table  4.  Comparison  of  the  Conditional  POF  Sensitivity  Equations 
with  Finite  Difference  Estimates 


Random  Variable  Parameter 

Sensitivity  Method 

Finite  Difference* 

Propagation  Scatter 

Median 

-5.148E-2 

-5.195E-2 

Propagation  Scatter 

COV 

4.765E-3 

5.500E-3 

Stress  Scatter 

Median 

2.037E+2 

2.093E+2 

Stress  Scatter 

COV 

9.042E-2 

9.860E-2 

Exceedance  Curve 

a 

min 

6.497E-3 

6.41 6E-3 

Exceedance  Curve 

a 

max 

3.260E-10 

3.602E-10 

*2E6  samples  used  for  all  analyses 
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Table  5.  Comparison  of  Disk  POF  Sensitivity  Equations 
with  Finite  Difference  Estimates 


Random  Variable  Parameter 

Sensitivity  Method 

Finite  Difference* 

Stress  Scatter 

Median 

1.035E+0 

1.066E+0 

Stress  Scatter 

COV 

7.409E-4 

7.796E-4 

Exceedance  Curve 

a 

min 

8.883E-9 

6.580E-8 

Exceedance  Curve 

a 

max 

5.989E-12 

5.258E-12 

Propagation  Scatter 
Median  -  cumulative 

-2.369E-4 

-2.553E-4 

Propagation  Scatter 

COV  -  cumulative 

5.891  E-5 

4.938E-5 

*  Importance  sampling 


4.4.4  Concluding  Remarks  on  Probabilistic  Sensitivities 

The  score  function  method  has  been  modified  to  fit  the  zone-based  methodology  in 
DARWIN.  The  score  function  method  has  advantages  in  that  the  sensitivities  are 
computed  without  any  additional  sampling  and,  therefore,  are  obtained  at  negligible 
additional  cost.  This  is  in  stark  contract  to  finite  differencing  which  requires  a  separate 
analysis  for  each  sensitivity.  Variance  estimates  have  been  obtained  such  that  the  lower 
and  upper  95%  confidence  bounds  for  each  sensitivity  is  output. 

Since  the  completion  of  the  probabilistic  sensitivity  information  into  DARWIN  under 
DUST  funding,  a  number  of  continuing  enhancements  have  occurred,  most  under  parallel 
USAF  funding.  These  enhancement  are  not  implemented  into  DARWIN  but  are  relevant 
for  fatigue  lifing  and  probabilistic  analysis  in  general. 

Using  the  same  score  function  concept  the  sensitivity  of  the  response  mean  and  response 
standard  deviation  can  be  computed  as 

%  =  £[Z(X)  Mi)]  (26) 

00 
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and 


^-  =  (£[Z(x)!k,(x)]-2//z£[Z(x)  k»(x)])/(2o-z)  (27) 

06 

where  Z  denotes  the  response,  e.g.,  cycles-to-failure. 

Also,  kernel  functions  for  correlated  normal  random  variables  have  been  developed. 


Figure  81.  Zone  Discretization  of  Rectangular  Cross  Section  with  Blow-up  of  Zone  2 
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4.5  Bayesian  Updating 


The  objective  for  this  task  was  to  demonstrate  the  fusion  of  data  from  modeling  (fatigue 
crack  growth  model)  and  testing  (NDE  estimate  of  crack  size). 

The  mathematics  for  combining  the  information  from  the  model  and  sensor  regarding  the 
crack  size  PDF  are  based  on  Bayesian  updating.  The  fundamental  equation  is 

f"(x)  =  kf(x)-f(x)  (28) 

where  / (x)  is  the  likelihood  PDF,  /'(x )  the  prior  PDF  ,  fix)  the  posterior  PDF,  and  k  is 
a  normalizing  constant  such  that  fix)  integrates  to  one,  i.e.,  k  =  \i\  fix)  ■  f  (x)dx . 


For  our  analysis,  /  (x)  represents  the  sensor  PDF,  fix)  represents  the  fracture  mechanics 
model-based  crack  size  PDF  after  a  flight  and  fix)  denotes  the  updated  PDF  that 
combines  the  information  from  the  model  and  sensor,  fix)  denotes  the  crack  size  that 
will  be  used  as  a  starting  crack  size  for  the  next  flight.  Equation  28  is  rewritten  specific  to 
our  analysis  as 


ffix)  =  kf+Gix)-fix)  (29) 

where  the  subscript  “iff  ’  represent  the  fracture  mechanics  model,  “M+G”  stands  for  the 
model  plus  crack  growth,  “S”  denotes  the  sensor,  and  i  denotes  the  flight  number.  After 
the  ith  flight,  the  crack  size  from  the  model  is  fM+Gix)  and  the  sensor  estimate  of  the  crack 
size  is  fix).  From  Equation  29,  the  starting  crack  size  for  the  simulation  of  the  next 
flight  using  the  fracture  mechanics  model  is  f'fix).  This  process  is  repeated  after  each 
flight. 

A  synopsis  of  the  methodology  is  shown  in  Figure  82;  the  crack  size  PDFs  are  fused 
through  a  Bayesian  updating  equation.  The  result  of  the  fusion  becomes  the  input 
distribution  to  the  model  for  the  next  flight. 
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crack  area  pdf  (mil2) 


Red  Curve  is  crack  size  PDF  predicted  from  the  fracture  mechanics  model, 
blue  curve  is  the  prediction  of  the  PDF  by  the  sensor,  green  is  the  resulting 
fused  PDF  the  becomes  and  input  to  the  next  flight. 


Figure  82.  Notional  Application  of  Bayesian  Updating 


4.5.1  Demonstration 

The  implementation  of  Bayesian  updating  was  accomplished  through  a  combination  of 
Matlab  (Bayesian  updating)  and  DARWIN  (crack  growth  analysis).  The  implementation 
is  nonparametric,  that  is,  no  assumptions  are  made  about  the  form  of  the  crack  size 
distributions  (either  DARWIN  or  the  sensor).  The  Monte  Carlo  crack  size  results  from 
DARWIN  were  output  to  a  file  and  analyzed  using  a  density  estimation  method  to  generate 
f‘M+G(x).  This  density  was  fused  with  the  sensor  density  fls(x)  to  produce  see 

Equation  29.  This  density  was  then  discretized  numerically  and  input  to  DARWIN  for  the 
next  flight. 

Some  general  notes  about  Bayesian  updating:  (1)  there  is  no  resulting  density  unless 
distributions  overlap,  and  (2)  the  density  with  the  tightest  variance  dominates.  That  is,  the 
“tightest”  density  contains  the  most  information  about  the  probable  location  of  the  crack 
size.  A  “non-informative”  prior  density  that  contains  no  information  on  the  probable  crack 
size  outside  of  the  given  bounds  would  be  a  uniform  distribution.  Therefore,  the  resulting 
distribution  after  Bayesian  updating  more  closely  takes  the  shape  of  the  tighter  distribution. 

The  sequential  analysis,  grow  the  defect  N  flights,  fuse  the  model  and  sensor  PDFs,  and 
input  the  fused  data  into  the  model  for  the  next  flight,  has  been  automated  through  a 
combination  of  Matlab  and  DARWIN,  see  Figure  83.  The  implementation  makes  no 
assumptions  about  the  DARWIN  analysis,  that  is,  whatever  content  is  contained  in  the 
DARWIN  template  file  will  get  executed.  This  provides  great  flexibility  regarding  the 
analyses  that  can  be  considered.  The  Bayesian  implementation  modifies  only  the  initial 


96 


crack  size  distribution  entered  into  DARWIN  and  reads  the  crack  size  results  after  a 
DARWIN  analysis. 


ec 


x 

o 

2 


Figure  83.  Schematic  of  Bayesian  Updating  Procedure 


4.5.2  Demonstration  Problem 

The  analysis  procedure  is  demonstrated  using  a  notched  beam  in  pure  bending  containing 
two  thin-fdm  magnetostrictive  ultrasonic  sensors.  The  crack  was  grown  in  fatigue  and 
interrogated  using  the  sensors. 

4.5.3  Initial  Crack  Size  Distribution 

A  central  issue  in  a  probabilistic  analysis  is  the  initial  crack  size  distribution.  Under  FAA- 
sponsored  efforts,  the  original  equipment  manufacturers  (OEMs)  have  developed  a  crack 
size  distribution  for  machined  holes  in  disks  which  represents  results  that  can  be  expected 
for  a  fleet  of  disks,  see  Figure  84  [42],  The  probability  density  is  “exponential-like”  in  that 
it  has  a  high  probability  of  small  cracks  with  a  long  right  tail  to  capture  the  possibility  of 
remote  large  cracks.  This  distribution  is  used  as  a  starting  distribution  for  the 
demonstration  problem.  As  will  be  shown,  the  initial  crack  size  distribution  is  modified 
significantly  based  upon  the  sensor  crack  size  estimation. 
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Figure  84.  Initial  Crack  Size  Distribution 


4.5.4  Thin-Film  Sensor  (TFS) 

The  signal-to-noise  ratio  from  the  TFS  was  compared  against  known  crack  area  obtained 
from  observation  of  post-inspection  marker  bands  on  the  specimen.  The  regression  results 
are  shown  in  Figure  85.  The  regression  results  indicate  a  constant  standard  deviation  of 
approximately  2349  mil2.  For  analysis  purposes,  the  TFS  result  was  truncated  at  a  crack 
size  of  2670  mil2  (S/N  =  1),  i.e.,  no  size  information  is  provided  below  this  crack  size. 

This  truncated  value  was  set  at  -2  standard  deviations  which  set  the  mean  size  as  7370  mil2 
(a  =  c  =  66  mil).  In  summary,  the  crack  size  density  estimate  from  the  sensor  was  modeled 
as  a  normal  distribution  with  mean  7370  mil2 ,  standard  deviation  2349  mil2  and  truncated 
in  the  left  tail  at  2670  mil2. 
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Wireless  -  With  Load 


Figure  85.  Thin-Film  Sensor  Experimental  Regression  Results 


4.5.5  Analysis  Results 

Constant  Sensor  Reading 

Various  scenarios  can  be  considered  using  the  developed  analysis  tools.  Figure  86  shows 
the  results  as  a  function  of  flights  with  the  sensor  density  unchanged  with  the  number  of 
flights.  The  results  show  the  model  PDF  becoming  more  narrow  and  approaching  the 
mode  of  the  sensor.  Most  interestingly  is  the  plot  summarizing  the  behavior  of  the  mean 
and  standard  deviation  of  the  fused  PDF.  The  results  from  flight  1  represent  the  initial 
crack  size  distribution  grown  1  flight.  After  fusing  with  the  sensor,  the  mean  value  jumps 
significantly  from  approximately  1 800  mil2  to  6000  mil2  because  the  sensor  has  found  an 
indication  of  a  defect.  Note  however,  that  the  standard  deviation  reduces  from  3500  mil2 
to  2000  mil2.  In  other  words,  we  have  a  much  stronger  belief  in  the  crack  size  after  the 
sensor  detects  an  indication. 

Sensor  Drift  to  the  Right 

Figure  87  shows  the  results  of  a  simulation  where  the  sensor  reading  is  taken  every  1000 
flights  and  the  sensor  drifts  to  the  right  by  400  mil2  per  every  1000  flights.  In  this 
scenario,  the  model  initially  under  predicts  the  PDF  relative  to  the  sensor  but  later  over 
predicts  relative  to  the  sensor.  The  mean  and  standard  deviation  show  a  similar  behavior 
as  discussed  previously;  the  mean  has  a  large  jump  initially  then  increases  more  slowly  and 
the  standard  deviation  drops  significantly  after  the  first  flight  then  reduces  more  slowly. 
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Figure  86.  Result  of  Bayesian  Updating  of  Model  and  Sensor  -  Fixed  Sensor 
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Figure  87.  Result  of  Bayesian  Updating  of  Model  and  Sensor  - 
Sensor  Drift  400  mil2/1000  Flights 
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4.5.6  Summary 

The  mathematics  necessary  to  combine  information  on  crack  size  from  models  and  sensors 
is  contained  in  Bayesian  updating.  This  technology  provides  a  straightforward  way  to 
update  the  crack  size  PDF  given  data  from  both  sources. 

The  Bayesian  updating  algorithm  was  implemented  in  software  using  DARWIN  for  the 
crack  growth  analysis  and  Matlab  for  the  Bayesian  updating.  The  methodology  makes  no 
assumptions  on  the  DARWIN  model  nor  any  assumptions  on  the  crack  size  distribution 
type. 

Two  demonstration  problems  were  exercised,  one  with  a  fixed  sensor  reading  (no  change 
with  flight  number)  and  one  with  a  sensor  that  drifted  400  mil2  with  each  1000  flights.  The 
problems  demonstrated  the  large  reduction  in  standard  deviation  of  crack  size  from  the 
assumed  initial  damage  PDF  achieved  by  having  a  sensor  reading. 


4.6  Benefits  of  Parallel  Processing 

The  objective  of  this  task  was  to  perform  an  introductory  look  into  parallel  processing 
methods  for  DARWIN  and  to  see  what,  if  any,  low  cost  items  could  be  handled  under 
DUST. 

Probabilistic  fracture  mechanics  is  inherently  computationally  time  consuming.  Although 
numerous  efficiency  enhancements  have  been  implemented  into  DARWIN  over  the  years, 
new  capabilities  have  also  been  added  which  require  increased  computational  time.  Most 
importantly,  analysis  requirements  for  military  missions  (many  time  points),  prognosis 
(cycle -by-cycle  integration),  residual  stress  analysis,  etc.  are  exceptionally  demanding. 
Table  6  below  shows  a  schematic  of  real  scenarios  and  the  corresponding  computational 
time  using  a  desktop  personal  computer  purchased  circa  2003.  Although  other  analyses 
will  have  different  execution  times,  it  is  clear  that  an  improvement  in  execution  time  is 
needed.  Fortunately,  parallel  processing  offers  an  avenue  to  address  this  issue. 
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Table  6.  Analysis  Scenarios  Driving  the  Need  for  Increased  Computational 
Efficiencies  in  Probabilistic  Fracture  Mechanics  Analyses 


Topic 

Reason  for  Computational  Expense 

Time  to  execute  - 
desktop  PC 

Simulation  of  onboard 

sensor 

Cycle-by-cycle  fatigue  crack  growth, 
single  inspection  scenario,  Monte 

Carlo  sampling 

24  hrs  with  5000  samples 
(more  samples  needed). 
Multiple  inspection 
scenarios  required 

Bayesian  updating  of 
engine  prognostics 
using  a  virtual  crack 
sensor  concept 

Cycle-by-cycle  fatigue  crack  growth, 
Bayesian  updating  of  crack  size, 
neural  network  training,  Monte  Carlo 
sampling 

Similar  but  slightly  longer 
to  “Simulation  of  onboard 
sensor”,  see  row  1 . 

Probabilistic  risk 
assessment  of  fatigue 
disk  life  subject  to 
inherent  defects 

Multiple  zone  analyses  required  for 
inherent  defects,  Monte  Carlo 
sampling 

50-100  hrs  for  2000 
samples.  More  samples 
required. 

Residual  stress  effects 
on  probabilistic  lifing  of 
engine  disk  materials 

Multiple  probabilistic  fatigue  analysis 
required  to  implement  GCE  algorithm 

20  hrs  with  2000  MC 
samples,  only  1  zone 

4.6.1  Parallel  processing  approaches 

There  are  several  levels  of  parallel  processing  typically  organized  by  “granularity”  of  the 
method.  OpenMP  is  well-suited  for  fine  grained  parallelism  (low  amounts  of 
communication  time),  Message-Passing  Interface  (MPI)  is  well-suited  for  coarse  grained, 
and  grid  computing  for  multiple  independent  analyses. 

4.6.2  Multi-Threading  (OpenMP) 

OpenMP  is  a  parallel  programming  model  for  shared  memory  computers  (multi-processor 
machines  where  all  processors  can  see  all  memory).  OpenMP  is  administered  by  an 
independent  oversite  board  (www.openmp.org).  As  such,  it  is  now  a  standard  for  Fortran 
and  C  and  supported  by  multiple  compiler  vendors.  The  OpenMP  web  site  also  contains 
tutorial  information  and  there  are  books  available  [43], 

OpenMP  is  a  multi-thread  approach  to  parallel  programming  -  a  master  execution  thread  is 
forked  into  any  number  of  threads  to  execute  a  user-defined  sequence  of  code,  see 
Figure  88.  The  forking  and  joining  of  threads  is  determined  by  the  user  using  compiler 
directives. 

One  of  the  advantages  of  OpenMP  is  that  an  executable  can  be  parallelized  “piece  by 
piece”  over  time.  That  is,  compiler  directives  can  be  inserted  in  the  source  code  to  multi- 
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thread  sections  of  code  as  time  permits.  This  is  in  contract  to  MPI  which  often  requires  a 
whole  scale  modification. 

OpenMP  compiler  directives  are  designed  such  that  compilers  that  do  not  support  OpenMP 
see  the  directives  as  comments.  Thus,  there  is  no  decrement  to  behavior  for  compilers  that 
do  not  support  OpenMP.  The  compiler  directives  are  documented  in  the  OpenMP 
specification,  tutorials  and  books. 

The  limitation  of  OpenMP  is  that  it  is  designed  for  shared  memory  machines  and  the 
number  of  processors  on  shared  memory  machines  is  limited;  although  this  is  changing. 

For  example,  a  top  end  desktop  machine  5  years  ago  contained  2  processors.  Now,  4 
processor  machines  are  routine  and  8  processor  machines  are  common.  This  trend  is  likely 
to  continue;  however,  shared  memory  machines  with  hundreds  or  thousands  of  processors 
are  not  expected  in  the  foreseeable  future.  This  fact  puts  a  limit  on  the  scalability  of 
parallel  algorithms  using  OpenMP. 


Figure  88.  Example  of  OpenMP  Fork- Join  Thread  Process 


4.6.3  Message-Passing  Interface 

MPI  is  the  most  general  method  for  parallel  processing.  MPI  requires  only  local  memory 
for  each  processor  and  a  means  of  communicating  between  processors  to  share  information 
[44],  This  model  matches  many  inexpensive  “clusters”  and  maps  to  the  Internet  in  general. 

MPI  is  a  standard  governed  by  an  oversite  board  (www.mpi-forum.org).  MPI  is  a  library 
not  a  compiler  enhancement,  hence  use  of  MPI  requires  installation  of  a  library  on  each 
system. 

The  essence  of  MPI  is  the  message;  the  sender  packs  data  into  a  message,  sends  the 
message  to  a  specific  user  or  users,  each  user  agrees  to  receive  the  message,  then  unpacks 
the  data. 
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A  key  element  that  often  defines  the  usefulness  of  MPI  is  the  computation  to 
communication  time  ratio.  Packing  data,  sending,  then  receiving  data  obviously  consumes 
central  processing  unit  (CPU)  and  wall  time.  A  good  algorithm  will  optimize  the  time 
spent  computing  versus  communicating. 

Because  MPI  is  well  suited  to  distributed  memory  machines  or  clusters,  theoretically,  it 
can  scale  to  any  number  of  processors,  depending  upon  the  algorithm.  Hence,  well 
designed  and  highly  parallel  algorithms  can  take  advantage  of  hundreds  or  thousands  of 
processors. 

4.6.4  Grid  Computing 

Grid  computing  is  the  coarsest  method  and  involves  decomposing  the  problem  into  a 
number  of  independent  analyses.  As  defined  by  Wikipedia:  “Distributed”  or  “grid” 
computing  in  general  is  a  special  type  of  parallel  computing  which  relies  on  complete 
computers  (with  onboard  CPU,  storage,  power  supply,  network  interface,  etc.)  connected 
to  a  network  (private,  public  or  the  Internet)  by  a  conventional  network  interface,  such  as 
Ethernet.  This  is  in  contrast  to  the  traditional  notion  of  a  supercomputer,  which  has  many 
processors  connected  by  a  local  high-speed  computer  “bus.”  The  advantage  of  grid 
computing  over  MPI  or  OpenMP  is  that  open-source  software  exists  that  can  be  easily 
installed  and  used  to  coordinate  the  independent  analyses. 

There  are  several  grid  computing  “engines”:  CONDOR  (http://www.cs.wisc.edu/condor/), 
Sun  Grid  Engine  (http://gridengine.sunsource.net/),  Globus  (www.globus.org),  IBM’s  Grid 
and  Grow  (http://www-03.ibm.com/grid/gridandgrow.shtml),  etc.  The  web  site 
http://www.gridcomputing.com/  offers  more  information  about  grid  computing. 

4.6.5  Application  of  Grid  Computing 

Under  FAA  funding,  DARWIN  was  modified  such  that  the  zone-based  analysis  would  fit 
within  the  grid  computing  paradigm  through  a  domain-decomposition  approach.  In  effect, 
each  zone  is  mapped  to  a  different  computer.  In  practice,  to  reduce  communication  time, 
multiple  zones  may  be  bundled  into  one  analysis  and  passed  to  a  computer. 

Under  DUST  funding,  the  multi-input  file  capability  was  exercised  using  the  grid 
computing  program  CONDOR  from  the  University  of  Wisconsin  using  a  bank  of  standard 
personal  computers  running  the  Windows  operating  system  located  at  UTSA.  The 
approach  and  results  are  contained  in  the  paper  “Application  of  Parallel  Processing  to 
Probabilistic  Fracture  Mechanics  Analysis  of  Gas  Turbine  Disks”  presented  at  the  2004 
AIAA  Structures,  Dynamics,  and  Materials  conference,  see  the  Appendix. 

The  results  revealed  the  effectiveness  of  the  approach;  an  efficiency  of  77%  was  obtained 
for  a  problem  consisting  of  6250  zones,  with  the  zones  divided  into  121  different  input 
files  and  run  over  18  computers. 
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4.6.6  Summary 

Parallel  processing  is  a  critical  technology  for  applying  DARWIN  to  analyses  of  military 
scenarios.  Parallel  processing  paradigms  of  MPI,  multi-threading  (OpenMP)  and  grid 
computing  were  briefly  reviewed.  Several  grid  computing  examples  were  developed  using 
a  parallel  network  of  personal  computers  running  Windows  2000  and  the  CONDOR  public 
domain  program.  The  approach  was  effective,  yielding  an  efficiency  77%  was  obtained 
for  a  problem  consisting  of  6250  zones.  The  efficiency  is  largely  controlled  by  the  ratio  of 
the  communication/computation  time.  This  ratio  can  be  approximately  controlled  by  the 
user  by  setting  the  number  of  zones  written  to  an  file.  Our  opinion  is  that  the  best  choice 
for  the  number  of  files  is  a  small  multiple,  say  3  or  4,  times  the  number  of  computers  in  the 
network.  This  allows  some  dynamic  load  sharing  without  excessive  communication  time. 
However,  this  approach  does  not  address  the  requirements  for  many  military  missions  that 
require  large  execution  times  for  a  single  zone. 

4.6.7  Future  Efforts 

The  grid  computing  approach  is  effective  for  multiple  zones  yet  has  some  significant 
deficiencies.  Foremost,  the  method  cannot  be  applied  to  single  zone  problems.  As  shown 
in  Table  6,  cycle-by-cycle  crack  growth  with  inspection  needed  for  prognosis  combined 
with  large  loading  histories  indicate  that  even  single  zone  problems  may  be  very 
computationally  demanding  problems  and  cannot  be  addressed  with  grid  computing.  As 
such,  it  is  proposed  that  a  significant  effort  be  expended  to  enhance  DARWIN  with 
OpenMP  and  or  MPI  capabilities.  The  logical  application  of  these  methodologies  would 
be  to  parallelize  the  Monte  Carlo  sampling  algorithm.  OpenMP  is  a  logical  first  step  but  it 
is  limited  to  shared  memory  computers.  However,  the  recent  progress  in  microprocessors 
is  towards  single  chips  with  many  CPU  cores.  Thus,  implementation  of  OpenMP  may  not 
be  as  limiting  as  at  first  glance.  MPI  is  the  most  general  approach  and  can  scale  to  any  size 
however  the  level  of  effort  to  implement  and  debug  is  substantially  higher.  Therefore,  it  is 
recommended  that  the  OpenMP  and  MPI  paradigms  be  explored  for  implementation  into 
future  versions  of  DARWIN. 


4.7  Combining  Classical  and  Physics-Based  Diagnostics  and  Prognostics 

Two  distinctly  different  approaches  to  diagnostics  and  prognostics  have  been  pursued,  and 
continue  to  evolve,  within  the  technical  community: 

1)  Classical  diagnostics  and  prognostics  -  these  methods  are  largely  data-driven 
and  employ  signal  processing  of  on-board  data  to  define  baseline  operating 
conditions,  deviations  from  these  baselines  (faults),  from  which  fault  limits  are 
established  and  used  to  trigger  inspection  and/or  maintenance  actions. 


106 


2)  Physics-based  diagnostics  and  prognostics  -  these  methods  are  based  on  an 
underlying  physical  model  for  material  damage  or  structural  degradation  -  for 
example,  the  deterministic  and  probabilistic  fatigue  crack  nucleation  and  growth 
models  described  in  the  preceding  section  of  this  report.  The  damage  rate  and  final 
instability  (failure)  of  the  material  is  independently  measured  and  expressed  in  terms 
of  a  physical  “driving  force”  such  that  the  result  is  a  material  property  (that  can  in 
some  instances  also  depend  on  environment  and  loading  rate).  These  material 
properties  are  then  combined  with  finite  element  stress  analysis  of  the  component  of 
interest  to  predict  the  nucleation  and  growth  of  damage,  as  well  as  final  failure  of  the 
component. 

Each  of  these  approaches  has  certain  advantages  and  disadvantages.  For  example,  the 
physics-based  methods  have  the  advantage  of  being  more  rigorous  and  potentially  more 
accurate  -  particularly  for  prognosis.  However,  the  disadvantage  of  these  methods  is  that 
they  can  be  computationally  intensive  for  certain  applications,  and  thus  may  not  be 
amenable  to  real-time  component  assessments;  they  also  require  a  detailed  understanding 
of  the  damage  or  degradation  process,  which  is  not  always  known.  In  contrast,  classical 
diagnostics  and  prognostics  approaches  are  generally  faster,  since  they  employ  more 
efficient  probabilistic  methods  based  on  artificial  intelligence,  neural  networks,  and 
Bayesian  statistics  that  do  not  require  detailed  knowledge  of  the  physics.  Consequently, 
they  are  applicable  to  a  wide  range  of  systems  (mechanical,  electrical,  fluid),  as  well  as  to 
higher  levels  within  the  system  -  for  example,  a  system  or  sub-system  that  may  contain 
many  components.  However,  this  lack  of  an  underlying  physical  model  means  that 
classical  methods  must  be  carefully  trained  (calibrated)  for  each  specific  application,  and 
depending  on  the  availability  of  empirical  data  their  selected  fault  limits  may  not  be 
optimum  or  universal.  One  consequence  of  the  fuzzy  nature  of  the  fault  limit  is  that 
classical  methods  are  generally  more  diagnostic  than  they  are  prognostic. 

In  view  of  the  above  trade  offs  between  the  two  approaches,  it  would  seem  that  some 
combination  of  these  methods  would  be  optimum  for  diagnostics  and  prognostics.  Thus, 
the  goal  of  this  task  was  to  assess  the  feasibility  and  potential  benefits  of  combining 
classical  diagnostic  and  prognostic  methods  with  physics-based  methods.  The  specific 
approach  in  this  case  was  to  look  for  turbine  engine  diagnostic  and  prognostic  problems 
that  might  demonstrate  enhanced  solutions  by  combining  the  capabilities  of  Smiths’ 
Probabilistic  Diagnostic  and  Prognostic  System  (ProDAPS)  and  SwRI’s  probabilistic 
fracture  mechanics  software  (DARWIN).  A  common  point  of  departure  for  this  exercise  is 
that  both  methods  recognize  the  considerable  uncertainty  that  can  exist  in  diagnostics  and 
prognostics,  thus  both  methods  are  probabilistically  based. 

ProDAPS  was  developed  on  a  parallel  DUST  Program.  ProDAPS  was  used  to  develop  and 
demonstrate  advanced  Prognostic  Health  Management  (PHM)  capabilities  for  a  range  of 
engine  and  aircraft  applications.  All  ProDAPS  elements  were  built  as  software 
components  and  were  designed  in  an  open  architecture  to  enable  technology  insertion  into 
both  airborne  and  ground-based  systems  for  existing  and  future  applications.  ProDAPS 
uses  advanced  artificial  intelligence  tools  to  provide  a  probabilistic,  knowledge  rich, 


107 


framework  to:  Integrate  all  available  types  of  health,  usage,  maintenance  and  reliability 
data;  Provide  advanced  data  analysis  capabilities  to  maximize  the  available  information; 
Fuse  and  reason  with  data  from  all  sources  to  provide  accurate  diagnostics  and  prognostics; 
Support  the  user  in  optimum  decision  making,  and;  Adapt  and  improve  its  performance 
through  the  mining  of  historical  data  to  discover  new  knowledge 

Several  approaches  to  assessing  the  effectiveness  of  ProDAPS  and  DARWIN  capabilities 
were  pursued.  The  goal  of  one  of  these  efforts  was  to  assess  if  ProDAPS  could  detect  fault 
events  that  would  be  potentially  detrimental  to  the  operating  environment  of  rotor/disks 
and  pass  them  to  DARWIN  for  more  effective  treatment  of  actual  engine  usage  and  its 
influence  on  component  fatigue  life.  The  output  from  DARWIN  could  also  be  fed  back  to 
ProDAPS  for  further  prognosis  and  decision  support.  This  exercise  was  focused  on 
detection  stalls  and  over-speeds  since  these  faults  should  be  detrimental  to  rotor/disk  health 
and  component  life.  This  exercise  required  that  data  be  fused  from  two  separate  data 
bases:  1)  the  Engine  Diagnostic  Unit  (EDU)  from  the  Comprehensive  Engine  Tracking  and 
Diagnostic  System  (CETADS),  which  among  other  things,  records  transient  data  when 
fault  events  on  engines  are  detected;  and  2)  data  from  Flight  Data  Recorders  (FDRs), 
which  among  other  things,  records  usage  in  terms  of  engine  RPM  vs.  time  histories  as  part 
of  the  Aircraft  Structural  Integrity  Program  (ASIP).  The  latter  data  are  required  to 
compute  fatigue  crack  growth  lives  using  DARWIN. 

Although  a  total  of  192  stall  events  and  1  over-speed  event  were  indicated  in  the  CETADS 
data  base,  we  were  able  to  find  the  corresponding  FDR  full  usage  data  for  one  stall  event. 
This  is  due  in  part  to  the  fact  that  FDRs  are  only  on  20%  of  the  fleet.  Moreover,  FDR  data 
sets  are  identified  by  aircraft  tail  number  (not  engine  serial  number),  and  CETADS  data  is 
only  identified  by  the  download  date,  which  is  not  necessarily  the  flight  date.  In  addition, 
further  scrutiny  of  the  data  determined  that  of  the  192  indicated  stall  events  only  one  was 
found  to  be  a  valid  stall,  and  there  was  no  matching  FDR  data  for  this  event.  Thus  no 
correlated  data  were  found  to  meet  the  exercise  objective.  The  major  conclusion  was  that 
extracting  correlated  data  from  CETADS  and  FDRs  is  both  inefficient  and  non-productive 
because  of  lack  of  common  identifiers  for  engines,  the  limited  number  of  aircrafts  that 
contain  FDRs,  as  well  as  the  inaccuracies  that  exist  in  the  EDU/CETADS  data. 

As  indicated  previously  in  Section  4.3,  statistical  analysis  of  usage  data  is  hampered  by  the 
fact  that  mission  type  is  often  not  recorder  in  the  FDR  database,  thereby  providing  the 
motivation  to  develop  an  automated  procedure  from  identifying  the  mission  given  the 
usage  (RPM  vs.  time).  Thus,  an  attempt  was  made  to  develop  a  data-driven  Mission 
Classifier  using  cluster  analysis.  Cluster  analysis  is  a  data  mining  tool  that  can  reveal 
associations  and  structure  in  data  that  is  not  initially  evident.  Using  statistical  analysis  this 
method  can  sort  cases  (e.g.  missions,  people,  things,  events)  into  groups,  or  clusters,  so  that 
the  degree  of  similarity  is  strong  between  group  members  and  weak  among  non-members. 
Cluster  analysis  was  assessed  for  its  ability  to  define  a  formal  classification  scheme  for 
mission  types  based  on  analysis  of  the  characteristics  of  flight  histories. 
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The  specific  data  employed  in  this  exercise  was  from  FI 00-229  engines  on  F-15  aircraft  at 
Lakenheath,  Nellis,  Mountain  Home,  and  Elmendorf  between  2002  and  2003.  Although 
FDR  data  was  available  from  nine  different  mission  types,  the  analysis  focused  on  the  two 
missions  -  Live  Fire  and  Instruments  and  Navigation.  The  data  metric  used  to  classify  the 
missions  was  the  percentage  of  flight  time  spent  at  various  Power  Lever  Angles  (PLAs)  as 
follows:  IDLE  to  18°,  18°-30°,  30°-50°,  50°-83°,  and  greater  than  83°.  Seven  different 
cluster  models  were  generated  having  two  and  eight  clusters  in  each  model.  Using  the 
Bayesian  Information  Criterion  (BIC),  the  best  cluster  model  for  both  mission  types  was 
determined  to  be  the  model  with  five  clusters.  However,  the  results  showed  no  clear 
distinction  between  the  to  different  mission  types.  Thus,  the  number  of  PLA  intervals  was 
reduced  from  seven  to  four  as  in  the  industry  definition  of  TACs.  In  this  case,  the  Live 
Fire  missions  were  best  described  by  a  cluster  model  with  five  clusters,  while  the 
instrument  and  navigation  missions  were  best  characterized  by  a  cluster  model  with  three 
clusters.  However,  as  shown  in  Figure  89,  it  was  found  that  the  characteristics  of  the 
above  cluster  models  were  very  close;  thus,  a  clear  characterization  parameter  to 
distinguish  between  the  two  mission  types  did  not  emerge. 

It  is  hypothesized  that  the  lack  of  success  in  using  the  above  cluster  analysis  to 
automatically  identify  missions  types  from  usage  is  due  either  to  the  relative  lack  of  data 
for  each  mission  type,  or  to  the  fact  that  the  missions  may  have  been  improperly  identified 
since  these  designations  are  currently  manually  logged  by  the  pilot  following  the  flight  -  a 
subjective  process  that  may  be  prone  to  error  or  inconsistencies. 


IDLE  CRUISE  INTERMEDIATE  MAX 


□  MIN_MIN 

■  MAXJVIAX 

■  Cluster  Model  5  Cluster  1 
Cluster  Model  5  Cluster  2 

■  Cluster  Model  5  Cluster  3 

■  Cluster  Model  5  Cluster  4 

■  Cluster  Model  5  Cluster  5 

■  MIN-MIN 

□  MAX-MAX 

■  cluster  model  3  cluster  1 
cluster  model  3  cluster  2 

□  cluster  model  3  cluster  3 


The  ordinate  in  the  figure  is  Bayesian  Information  Criterion  (BIC)  and  the  abscissa  shows 
intervals  of  Power  Lever  Angle  (PLA)  per  the  industry  definition  of  Total  Accumulated  Cycles. 


Figure  89.  Cluster  Analysis  Indicated  That  Cluster  Model  Five  Best  Characterized  the 
Live  Fire  Missions,  and  Cluster  Model  3  Best  Characterized  the  Instrument  and 
Navigation  Missions;  However,  These  Characteristics  Were  Not  Sufficiently 
Different  to  Uniquely  Discriminate  Between  the  Two  Mission  Types 


109 


The  above  analyses  served  to  highlight  two  important  issues:  1)  the  limitations  in  using 
current  data  bases  for  the  tasks  performed  here  -  specifically,  data  are  not  consistently 
cross  referenced  from  one  data  base  to  the  next  so  that  they  can  be  fused  for  a  given 
mission,  and  2)  the  difficulties  encountered  in  attempting  to  develop  an  automated  mission 
classifier  stimulated  thought  that  led  to  the  eventual  development  of  the  Probabilistic 
Mission  Identification  method  described  in  Section  4.3. 

Regarding  the  need  for  data  strategies  that  facilitate  sharing  and  processing  data  across 
data  bases  and  organizational  boundaries,  the  current  situation  not  only  limited  the  above 
tasks,  but  would  also  limit  diagnostics  and  prognostics  in  general,  as  well  as  the  broader 
Department  of  Defense  (DoD)  goals  under  the  Condition  Based  Maintenance  Plus  (CBM+) 
initiative  [45].  Fortunately,  this  limitation  has  been  recognized  in  the  guidelines  that  are 
being  prepared  to  implement  these  important  initiatives  [46], 
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5.  Results  and  Discussion  on  Technology  Transition:  New  DARWIN 
Capabilities 

5.1  Summary  of  DARWIN  Enhancements 

The  following  new  DARWIN  features  were  implemented  under  this  project: 

•  FaNG  Model  (General  Inherent  Analysis  Mode,  ID  Surface  Damage  Analysis  Mode) 

•  Residual  Stress  Capability  (General  Inherent  Analysis  Mode) 

•  Mission  Mixing  (All  Analysis  Modes) 

•  Probabilistic  Sensitivities  (General  Inherent  Analysis  Mode). 

A  brief  description  of  the  implementation  and  use  of  these  features  is  provided  in  the 
sections  that  follow. 

5.2  Implementation  in  DARWIN 
5.2.1  FaNG  Model  Implementation 

A  new  DARWIN  capability  was  implemented  that  allows  the  User  to  select  the  FaNG 
model  for  crack  nucleation  and  growth  computations.  The  FaNG  model  predicts  the  crack 
nucleation  and  growth  lives  for  a  specified  initial  crack  size.  It  consists  of  (1)  a  smooth 
specimen  total  life  computation,  (2)  a  smooth  specimen  crack  growth  life  computation,  and 
(3)  a  surface  crack  in  plate  crack  growth  life  computation.  The  smooth  specimen  total  life 
calculation  is  performed  using  equations  that  are  defined  in  a  separate  module  (the 
DARWIN  “formation  module”).  This  module  is  called  during  execution  of  DARWIN,  and 
can  be  edited  by  the  User  as  necessary  for  definition  of  the  equations  associated  with 
smooth  specimen  total  life.  The  smooth  specimen  crack  growth  life  is  performed  within 
DARWIN  using  a  surface  crack  in  round  bar  stress  intensity  solution  (SC07)  that  was 
recently  implemented  specifically  for  this  purpose.  The  surface  crack  in  a  plate  crack 
growth  life  computation  is  performed  using  existing  stress  intensity  solution  SC  17.  Note 
that  the  FaNG  model  is  currently  limited  to  surface  crack  in  plate  geometries. 

During  execution,  the  smooth  specimen  nucleation  life  is  obtained  by  subtracting  the 
smooth  specimen  crack  growth  life  from  the  smooth  specimen  total  life.  For  smooth 
specimen  crack  nucleation  life  computations,  the  stress  value  is  set  to  the  value  at  the  User- 
specified  nucleation  (initial  crack)  depth.  This  stress  value  is  also  used  for  the  smooth 
specimen  crack  growth  life  computations,  where  a  uniform  stress  field  is  assumed.  A 
small  crack  correction  factor  is  applied  to  all  crack  growth  computations  to  account  for 
changes  in  the  crack  growth  rate  associated  with  small  fatigue  cracks. 
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5.2.2  Residual  Stress  Capability  Implementation 

A  new  DARWIN  capability  was  implemented  for  treatment  of  surface  residual  stress 
values.  This  new  capability  allows  the  User  to  define  a  residual  stress  field  that  appears  as 
a  rectangular  region  that  is  superimposed  on  the  2-dimensional  axisymmetric  finite 
element  model  that  appears  in  the  DARWIN  GUI.  Stress  values  within  the  region  that  are 
applied  to  crack  nucleation  and  growth  life  computations  are  based  on  the  sum  of  the 
residual  stress  values  and  the  finite  element  stress  values. 

5.2.3  Mission  Mixing  Implementation 

The  existing  DARWIN  mission  mixing  feature  provides  the  capability  to  assess  the  life 
and  probability  of  fracture  of  a  component  that  has  the  same  mission  mix  over  the  entire 
design  life.  For  example,  consider  a  mission  mix  that  consists  of  three  missions  (Missions 
A,  B,  and  C)  that  have  equal  mix  percentages.  If  the  existing  DARWIN  mission  mixing 
capability  is  used,  the  mission  sequence  “ABCABCABC  ...”  is  applied  repeatedly  over  the 
lifetime  of  the  component. 

Under  this  project,  the  DARWIN  mission  mixing  capability  was  enhanced  to  allow  the 
User  to  specify  changes  in  the  mission  mix  at  specified  flight  intervals.  For  example, 
consider  two  mission  mixes  (Mission  Mixes  1  and  2)  that  begin  at  different  flight  intervals. 
Mission  Mix  1  consists  of  missions  A,  B,  and  C,  and  Mission  Mix  2  consists  of  missions  D 
and  E.  Mission  Mix  1  begins  at  0  flight  cycles,  and  Mission  Mix  2  begins  at  5001  cycles. 
For  this  scenario,  the  mission  sequence  “ABCABCABC  . . .”  is  applied  repeatedly  over  the 
first  5000  cycles,  and  then  the  mission  sequence  “DEDEDE  ...”  is  applied  over  the 
remainder  of  the  life. 

5.2.4  Probabilistic  Sensitivities  Implementation 

The  sensitivity  methods  developed  by  Karamchandani  [47]  and  Wu  [48]  were  extended  for 
application  to  additional  parameters  of  random  variables.  Analytical  expressions  for 
probabilistic  sensitivities  were  developed  for  the  stress  scatter,  life  scatter,  and  initial 
anomaly  area  random  variables.  Variance  and  associated  confidence  bound  estimates  were 
also  developed  for  these  variables.  Probabilistic  sensitivities  and  confidence  bound 
expressions  were  implemented  in  DARWIN  for  these  three  random  variables.  A  significant 
feature  of  the  implementation  is  that  sampling-based  sensitivities  are  computed  directly 
from  the  Monte  Carlo  results  associated  with  failure  probability  computations,  and  are 
therefore  very  inexpensive  to  compute.  Complete  details  regarding  the  equations  are 
provided  in  [40], 
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5.3  New  DARWIN  GUI  Features 

5.3.1  Enabling  New  Features  in  DARWIN  GUI 


The  DARWIN  enhancements  associated  with  this  project  were  implemented  as  password 
protected  features.  To  enable  these  features,  the  User  must  provide  a  unique  password 
(Figure  90). 


1  ii  DARWIN  6.1.0  (Build  1295) 

File  View 

Tools 

Projects  Help 

i - 

a  a 


Set  license  key. 


<9  Convert  result  files... 

Edit  results  file... 
Recover  results  file. . . 


<9  Batch  report  generation. . . 
<9  Convert  geometry... 

Batch  processing... 

Create  parallel  file  set... 


Create  locked  library  files... 

<9  CAAM/Fleet  assessment  analysis  tool... 
Enable  optional  feature. 


Figure  90.  The  DARWIN  Enhancements  Associated  with  This  Project 
Were  Implemented  as  Password-Protected  Features 


113 


5.3.2  Use  of  FaNG  Model  in  DARWIN  GUI 


As  shown  in  Figure  91,  the  FaNG  model  was  implemented  as  a  life  option  associated  with 
the  DARWIN  Formation  module.  When  the  FaNG  option  is  selected,  the  User  must 
specify  values  for  the  ao  constant  (small  crack  correction)  and  the  diameter  of  the  specimen 
associated  with  smooth  specimen  crack  growth  life  (SC07  diameter). 


Figure  91.  The  FaNG  Model  Was  Implemented  as  a  Life  Option  Associated 
with  the  DARWIN  Formation  Module 


5.3.3  Use  of  Residual  Stress  Capability  in  DARWIN  GUI 


A  capability  for  treatment  of  residual  stress  (RS)  values  was  implemented  in  DARWIN. 
This  feature  is  activated  by  selecting  use  of  “residual  stress  profiles”  located  in  the 
DARWIN  “Setup”  menu  (Figure  92).  Definition  of  residual  stress  profiles  in  DARWIN 
consists  of  the  following  steps: 


(1)  Create  an  RS  region.  Click  the  _  button  in  the  Geometry  menu  to  bring  up  the 


Stress  Region  Editor  window,  (Figure  93).  In  the  Stress  Region  Editor  window, 

click  the  1  ¥  I  button  to  create  a  stress  region.  The  RS  region  is  defined  by  a  region 
edge  (width  of  the  RS  region)  and  values  of  the  RS  versus  physical  distance  from 
the  starting  edge  (length  of  the  RS  region). 

(2)  Draw  Region  Edge.  Click  the  |  ND,aw rea'°n  ed°e  |  button  to  draw  the  region  edge.  On 
the  finite  element  geometry,  click  at  the  starting  point  of  the  region  edge,  hold  the 
mouse,  drag  it  to  the  end  point  of  the  region  edge,  then  release  the  mouse. 
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Global  Factors 
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This  feature  is  activated  by  selecting  use  of  “Residual  Stress  Profiles”  located  in  the  DARWIN 
“Setup”  menu. 


Figure  92.  A  Capability  for  Treatment  of  Residual  Stress  Values 
Was  Implemented  in  DARWIN 


Figure  93.  The  DARWIN  Stress  Region  Editor  Allows  the 
User  to  Define  Residual  Stress  Gradient  Values 
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(3)  Edit  Region  Factors.  Enter  values  of  RS  vs.  distance  in  the  Region  Factors  table. 

Click  the  ®  button  to  insert  a  new  row  and  click  the  ®  button  to  delete  a  row. 

The  values  can  also  be  imported  from  a  text  file.  Click  the  button  to  open  the 

text  file  containing  region  factors.  To  export  region  factors,  click  the  IMI  button  to 
specify  the  filename.  The  region  factor  vs.  distance  plot  can  be  viewed  by  clicking 

the  ^  button. 


The  resulting  residual  stress  field  appears  as  a  rectangular  region  that  is  based  on  the  length 
of  the  residual  stress  gradient  and  the  length  of  the  user-specified  line  near  the  surface  of 
the  model.  The  region  is  superimposed  on  the  2-dimensional  axisymmetric  finite  element 
model,  and  stress  values  within  the  region  are  based  on  the  sum  of  the  residual  stress  and 
the  finite  element  stress  values  (Figure  94). 


5.3.4  Use  of  Mission  Mixing  in  DARWIN  GUI 

The  existing  DARWIN  mission  mixing  capability  was  enhanced  under  this  project  to  allow 
the  User  to  specify  changes  to  the  mission  mix  as  a  function  of  flight  cycles.  As  shown  in 
Figure  95,  a  mission  mix  scenario  menu  was  added  for  this  purpose.  The  capability  can  be 
used  in  conjunction  with  the  FaNG  and  residual  stress  models. 

To  use  the  enhanced  mission  mix  capability, 


1. 

2. 


3. 


Select  the  Mission  editor  tab  on  the  Missions  menu  to  define  missions. 

Select  the  Mission  mixer  tab  to  define  mission  mixes.  By  default,  DARWIN  creates  a 
mission  mix  (Mix  1)  with  the  first  mission.  Edit  the  mission  mix  by  adding/deleting 

missions  or  changing  the  order  of  missions.  Click  the  *  button  to  create  and  edit 
more  mission  mixes. 

Specify  the  mission  mixes  in  the  Mix  scenario  table  (located  at  the  bottom  of  the 
Mission  Mixer  tab).  Click  the 


<=[]□  Add  mix 


button  to  add  mission  mixes  to  the  table  and 


specify  the  starting  flight  cycle  for  each  active  mission  mix. 
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The  resulting  residual  stress  field  appears  as  a  rectangular  region  that  is  superimposed  on  the 
2-dimensional  axisymmetric  finite  element  model. 


Figure  94.  Definition  of  Residual  Stress  Profiles  in  DARWIN 
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Figure  95.  A  New  DARWIN  Capability  Was  Developed  That  Allows  the 
User  to  Specify  Changes  to  the  Mission  Mix 
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6.  Results  and  Discussion  on  Probabilistic  Simulations 

6.1  Mission  Mixing 

In  this  section,  the  DARWIN  mission  mixing  capability  is  described  briefly  and  illustrated 
for  representative  engine  usage  histories.  As  shown  in  Figure  96,  the  load  steps  associated 
with  each  mission  are  specified  using  the  DARWIN  Mission  Editor  menu.  The  mission 
mix  is  defined  by  specifying  the  sequence  and  number  of  repetitions  of  each  mission 
(Figure  97).  The  load  steps  associated  with  a  representative  mission  mix  are  shown  in 
Figure  98. 

The  DARWIN  mission  mixing  capability  was  illustrated  for  the  representative  gas  turbine 
engine  hub  component  shown  in  Figure  99.  Raw  RPM  values  were  obtained  from  engine 
flight  data  recorder  data  for  several  F-16/F100  missions.  The  RPM  values  were  converted 
to  stress  values  using  an  empirical  algorithm  and  sorted  into  min-max  pairs  using  the 
DARWIN  rainflow  cycle  counting  algorithm. 


Normalized  probability  of  failure  values  associated  with  representative  engine  usages  are 
shown  in  Figure  100.  In  Figure  101,  the  influence  of  mission  mix  ratio  on  normalized 
probability  of  failure  values  is  shown  for  the  representative  mixes  indicated  in  Table  7.  It 
is  interesting  to  note  that  when  the  mix  percentages  are  equal  (Mix  1),  the  probability  of 
failure  is  not  equal  to  the  average  of  the  probability  of  failure  values  of  the  individual 
missions.  This  is  probably  due  to  the  increased  K  values  associated  with  large  crack  sizes 
associated  with  the  most  damage  mission  (similar  to  sequence  effects).  As  expected,  the 
probability  of  failure  is  greater  for  mission  mixes  with  a  higher  percentage  of  damaging 
missions  (Mixes  2  and  3). 
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Figure  96.  Specifying  the  Load  Steps  Associated  with  a  Mission  in  DARWIN 
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Figure  97.  Specifying  the  Mission  Mix  in  DARWIN 


Mssion  1  12  2  1 


Figure  98.  A  Representative  Mission  Mix 
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Figure  99.  The  DARWIN  Mission  Mixing  Capability  Was  Illustrated 
for  a  Representative  Gas  Turbine  Engine  Hub  Component 


Table  7.  Mission  Mix  Percentages  for  Illustrative  Example 


Mission  ID 

Mix  1 

Mix  2 

Mix  3 

Live  Fire 

25% 

50% 

16.7% 

Instruments  &  Navigation 

25% 

16.7% 

50% 

ACBT 

25% 

16.7% 

16.7% 

Live  Fire  Missiles 

25% 

16.7% 

16.7% 
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Flights  (Thousands) 


Figure  100.  Normalized  Probability  of  Failure  Values  Associated 
with  Representative  Engine  Usages 


Figure  101.  Influence  of  Mission  Mix  Percentage  on  Normalized  Probability 
of  Failure  Values  for  Several  Representative  Mixes 
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6.2  Forecasting  Reliability 


Previous  versions  of  DARWIN  provided  a  framework  for  forecasting  the  reliability  of  a 
component  that  has  the  same  mission  mix  over  the  entire  design  life.  As  described  in 
Section  5,  DARWIN  was  enhanced  under  this  project  to  allow  the  User  to  specify  changes 
in  the  mission  mix  at  specified  flight  intervals.  This  enhancement  extends  the  reliability 
forecasting  capability  for  potential  application  to  fleet  management  decisions,  because  it 
allows  the  analyst  to  predict  the  influence  of  changes  in  applied  usage  on  the  overall  risk  of 
fracture. 

The  DARWIN  mission  mixing  enhancement  was  illustrated  for  the  representative  gas 
turbine  engine  hub  component  shown  in  Figure  99  (see  Section  6.1  for  further  details). 
Three  mission  mix  scenarios  were  considered  as  indicated  in  Table  8.  Although  DARWIN 
can  be  used  for  assessment  of  large  numbers  of  missions  and  mission  scenarios,  these 
simple  scenarios  were  selected  to  demonstrate  the  maximum  influence  of  changes  in  a 
single  mission  type  over  the  service  life  of  a  component.  The  results  are  shown  in 
Figure  102.  For  this  example,  it  can  be  observed  that  the  probability  of  fracture  is  usually 
largest  for  mix  scenario  2,  in  which  the  most  damaging  mission  mix  is  applied  during  the 
first  half  of  the  design  life.  Note  that  DARWIN  currently  includes  only  a  single  stress 
scatter  random  variable,  so  the  stress  scatter  values  associated  with  the  mix  scenarios  were 
based  on  an  average  of  the  values  associated  with  the  individual  missions  within  each  mix. 


Table  8.  Mission  Mix  Scenarios  for  Illustrative  Example 


Mix  Scenario 

Mix  1 

0  to  4000  flights 

Mix  2 

4000  to  8000  flights 

i 

50%  Live  Fire 

50%  Live  Fire 

50%  ACBT 

50%  ACBT 

2 

100%  Live  Fire 

100%  ACBT 

3 

100%  ACBT 

1 00%  Live  Fire 
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Flights  (Thousands) 


Figure  102.  Influence  of  Mission  Mix  Scenario  on  Normalized  Probability  of 
Failure  Values  for  Several  Representative  Mix  Scenarios 


6.3  Probabilistic  Simulation  of  On-board  Monitoring 

The  objective  of  this  task  was  to  compare  the  effectiveness  of  an  accurate  midlife  depot 
inspection  versus  a  coarser  continual  (once  per  flight)  onboard  inspection. 

The  traditional  approach  to  inspection  simulation  is  to  assume  each  inspection  is 
independent  from  the  previous.  This  assumption  is  appropriate  for  the  traditional 
application  of  inspections,  i.e.,  different  operators,  equipment,  and  location,  and  widely 
separated  by  the  time  of  inspection.  Well-established  computer  models  for  determining 
probability-of- fracture  (POF)  estimates  for  aircraft  structures  [49];  and  gas  turbine  rotor 
disks  [50]  consider  infrequent  or  opportunity  inspections,  and;  therefore,  model  the 
inspections  as  independent. 

However,  the  assumption  of  independent  inspections  leads  to  dramatic  predictions  of  non- 
realizable  efficacy  for  continual  inspections  and  the  results  are  frequency  dependent; 
almost  any  size  crack  can  be  detected  by  inspecting  more  and  more  frequently,  even 
without  any  crack  growth.  The  assumption  of  independent  inspections  is  severely  non¬ 
conservative  and  provides  a  false  estimate  of  the  accuracy  of  continual  inspections. 

Conversely,  it  is  asserted  and  discussed  below  that  continual  inspections  should  not  be 
treated  as  independent  from  previous  inspections  but  as  dependent;  since  the  same  sensor 
is  performing  subsequent  inspections  on  the  same  defect,  the  sensitivity  of  the  sensor  in 
detecting  the  crack  should  be  the  same  as  for  previous  inspections,  or  nearly  so. 
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Dependant  inspections  enforce  that  the  sensitivity  of  a  sensor  at  one  inspection  is  the  same 
at  subsequent  inspections. 

An  examination  of  redundant  fluorescent  penetrant  inspection  (FPI)  substantiated  the  fact 
that  repeated  inspections  should  not  be  treated  as  independent  and  that  doing  so  is 
nonconservative  [51].  In  this  study,  specimens  were  inspected  using  the  same  FPI  method 
by  two  inspectors  of  approximately  the  same  skill  level.  The  resultant  probability-of- 
detection  from  both  inspections  was  quantified  using  data  and  compared  with  the 
independent  assumption.  The  data-based  Probability  of  Detection  (POD)  considering  both 
inspections  is  significantly  lower  than  the  independent-based  POD  (which  can  be  obtained 
by  a  multiplication  of  one  minus  the  probabilities).  As  stated  [51],  “The  fluorescent 
penetrant  inspection  process  is  not  independent  inspection-to-inspection  and  therefore  the 
probability  of  detection  for  redundant  FPI  cannot  be  obtained  by  a  simple  multiplication  of 
probabilities.  The  correct  POD  in  a  multiple  inspection  system  is  determined  assuming 
dependence.” 

Mathematical  bounds  prove  that  independent  and  dependent  inspections  provide  the  upper 
and  lower  bounds,  respectively,  on  the  cumulative  probability  of  detecting  a  defect.  That 
is,  assuming  dependent  inspections  is  a  conservative  assumption.  Intermediate  results 
between  the  independent  and  dependent  assumptions  can  be  obtained  by  correlating  the 
inspections  with  a  correlation  coefficient  between  zero  (independent)  and  one  (dependent). 
The  sum-of-uniforms  method  [52,  53]  is  ideal  to  implement  any  amount  of  correlation 
between  inspections. 

Given  the  assertion  of  dependent  inspections,  it  is  shown  that  the  efficacy  of  an  on-board 
sensor  is  principally  determined  by  the  probability  of  detecting  a  crack  of  critical  size  and, 
by  and  large,  unaffected  by  the  shape  of  the  POD  curve.  This  information  can  be  used  to: 
a)  determine  the  effectiveness  of  implementing  a  particular  sensor  on  a  structural 
component,  or  b)  determine  the  necessary  design  specification  a  sensor  must  satisfy  to 
make  an  implementation  effective  for  component  safety  and  management  of  the  asset. 

6.3.1  Numerical  Example 

The  problem  consisted  of  a  surface  crack  on  the  bore  of  a  gas  turbine  disk.  The  stress 
spectra  applied  to  the  crack  is  derived  from  RPM  values  versus  time  obtained  from  a  flight 
data  recorder  representative  of  an  Air-to-Ground  military  mission.  The  mission  contains 
115  min-max  fatigue  loads  which  were  applied  per  each  flight  simulation.  The  initial 
crack  size  was  a  semi-circular  defect  with  an  initial  crack  depth  specified  by  an  exponential 
probability  distribution  and  a  1  -to- 1  aspect  ratio  with  a  mean  0.432  mm  (17  mils).  The 
initial  crack  size  distribution  was  representative  of  the  cracks  initiated  on  holes  in 
compressor  disks  due  to  machining.  Other  initial  defect  size  distributions  can  be 
considered  without  altering  the  main  conclusions.  Material  properties  for  a  titanium 
material  were  used.  The  material  fractured  at  a  critical  crack  surface  length  of 
approximately  15.24  mm  (600  mils)  under  the  applied  loading.  Table  9  provides  a 
summary  of  the  problem  inputs.  The  paper  in  the  Appendix  has  more  details. 
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Table  9.  Example  Problem  Definition 


Property 

Description 

Comment 

Loading 

Spectrum  loading 
representative  of  an  Air-to- 
Ground  mission 

11600  time  points,  231  load  pairs 

da/dN-AK 

C=5.248E-1 1  (9.25E-13  SI), 
m=3.87,  Kc=58.7  ksi-inA-1/2 
(64.6  MPa-mA-1/2) 

Paris  law  with  no  threshold  (C  -  constant,  m 
-  exponent,  Kc  -  fracture  toughness) 

US  -  units:  da/dN  (in/cycle),  AK  (ksi-inA-1/2) 

SI  -  units:  da/dN  (m/cycle),  AK  (MPa-mA- 
1/2) 

Initial  crack  size 

Probability  distribution 
representative  of  crack  at  holes 
of  compressor  disks 

Surface  crack  in  the  bore  of  a  gas  turbine 
disk.  1-1  aspect  ratio 

The  POD  curve  was  modeled  with  a  lognormal  distribution, 

POD(A  <  a)  =  0[ln(a/l)/cr]  (30) 

where  a  represents  crack  size,  X  is  the  median,  COV  is  the  coefficient  of  variation, 
c =  ln[l  +  COF2],  and  O  represents  the  standard  normal  cumulative  distribution  function. 

The  COV  used  for  the  POD  curves  is  chosen  as  thirty-two  percent  based  on  analysis  of 
existing  industry  ultrasonic  POD  curves. 

The  reduction  in  POF  obtained  by  a  single  but  precise  mid-life  inspection  was  compared  to 
the  results  from  less-precise  but  frequent  continual  inspections.  All  continual  inspections 
were  modeled  as  dependent;  therefore,  these  results  are  conservative.  Figure  103  shows 
the  normalized  POF  versus  flight  cycles  (the  results  are  normalized  to  the  POF  without 
inspection  at  8000  flight  cycles)  considering:  a)  no  inspection,  b)  a  single  mid-life 
inspection  with  a  median  POD  value  of  0.762  mm  (30  mils)  and  a  COV  of  32%,  c) 
continual  inspections  (once  per  flight)  with  a  median  POD  value  of  5.080  mm  (200  mils) 
and  a  COV  of  32%,  d)  continual  inspections  with  a  median  POD  value  of  10.16  mm 
(400  mils)  and  a  COV  of  32%,  e)  continual  inspections  with  a  median  POD  value  of 
15.24  mm  (600  mils)  and  a  COV  of  32%. 
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Figure  103.  Comparison  of  Continual  and  Single  Mid-Life  Inspections 


Figure  103  demonstrates  that  a  single  mid-life  inspection  reduces  the  POF  to 
approximately  22%  relative  to  the  POF  without  inspection  and  effectively  detects  almost 
all  defects  that  could  case  fracture  within  the  ensuing  4000  flight  cycles.  Continual 
inspections  with  a  sensor  that  exhibits  a  POD  of  median  value  equal  to  5.080  mm  (200 
mils)  and  a  COV  equal  to  32%  is  far  superior  than  a  single  more  precise  mid-life 
inspection,  reducing  the  POF  to  near  zero  even  though  the  precision  of  the  on-board  sensor 
is  far  less  than  the  depot  inspection.  Continual  inspections  with  a  sensor  that  exhibits  a 
POD  of  median  value  equal  to  10.16  mm  (400  mils)  and  a  COV  equal  to  32%  reduces  the 
POF  to  11%  relative  to  the  POF  without  inspection.  Continual  inspection  with  this  POD  is 
still  superior  to  a  single  more  precise  mid-life  inspection  through  8000  cycles  although  the 
slopes  indicate  that  at  some  number  of  flight  cycles  the  single  mid-life  inspection  will  be 
superior.  Continual  inspections  with  a  sensor  that  exhibits  a  POD  of  median  value  equal  to 
15.24  mm  (600  mils)  and  a  COV  equal  to  32%  reduces  the  POF  to  50%  relative  to  the  POF 
without  inspection  and  clearly  is  inferior  to  the  single  mid-life  inspection.  We  note  that  the 
critical  crack  size  for  this  problem  is  approximately  15.24  mm  (600  mils). 

For  dependent  inspections,  the  probability  of  detecting  a  flaw  before  fracture  is 
approximately  equal  to  the  value  of  the  POD  curve  at  the  critical  crack  size.  Therefore,  the 
reduction  in  POF  is  PF  reduced  ~  [1  -  POD(acnl)\  and  the  POF  with  inspection  can  be  closely 
approximated  given  the  POF  without  inspection  for  any  number  of  cycles  M  as 

Pf.«W)  *  PF,red»ced  '  PF^iM)  =  [1  -  POD^J]  ■  P^JM)  (31) 
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Equation  3 1  has  been  verified  by  considering  four  vastly  different  contrived  POD  curves 
developed  such  that  POD(acrit )  =  70%  for  each  curve,  see  Figure  104.  In  this  case 
Pf, reduced  =  1 _  POD(acrit)  ~  30%  for  each  POD  curve.  As  a  result,  the  POF  with  inspection 
was  approximately  the  same  regardless  of  which  POD  curve  was  used,  see  Figure  105. 
These  results  indicate  that  an  estimate  of  the  reduction  in  POF  due  to  continual  inspections 
can  be  ascertained  by  knowing  the  sensor’s  POD  value  at  the  critical  defect  size. 
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Figure  104.  Contrived  POD  Curves 
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Figure  105.  Reduced  POF  Resulting  from  Contrived  POD  Curves 


6.3.2  Summary 

During  this  research,  a  methodology  was  developed  and  implemented  to  compute  the 
reduction  in  POF  due  to  continual  (once  per  flight)  inspections  based  upon  the  concept  of 
on-board  crack  sensors.  This  issue  was  investigated  through  the  analysis  of  fracture  of  a 
gas  turbine  disk.  A  critical  outcome  of  the  research  was  that  accurate  simulation  of 
continual  inspections  requires  that  inspections  be  treated  as  dependent  rather  than  the 
traditional  modeling  approach  of  independent  inspection.  Dependent  simulation  is 
conservative  whereas  independent  simulation  may  be  grossly  unconservative. 

A  second  critical  outcome  is  that  the  reduction  in  POF  obtained  by  continual  inspections 
depends  almost  solely  upon  the  probability  of  detection  of  the  sensor  at  the  critical  crack 
size,  POD(acrit)-  This  information  can  be  used  to  determine  the  required  accuracy  for  an 
on-board  sensor  to  prove  effective.  This  assertion  was  verified  from  the  analyses  of  four 
drastically  different  POD  curves  with  the  same  POD  value  at  the  critical  crack  size;  the 
resulting  reduced  POF  from  the  analyses  were  all  the  same. 
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Finally  an  analysis  comparing  the  reduction  in  POF  from  a  precise  single  depot  inspection 
versus  various  much  coarser  but  continual  inspections  demonstrated  that  continual 
inspection  may  provide  a  significant  reduction  in  POF  compared  to  the  depot  inspection, 
even  if  the  onboard  sensor  is  several  times  less  sensitive  than  the  depot  inspection.  This  is 
due  to  the  fact  that  with  a  large  number  of  on-board  inspections,  multiple  opportunities 
exist  to  detect  larger  (but  non-critical  defects)  than  the  single  mid-life  inspection. 


6.4  Application  of  Bayesian  Updating  to  Health  Monitoring  and  Prognosis 

The  utility  of  the  Bayesian  updating  method  described  in  Section  4.5  was  demonstrated  by 
employing  it  to  fuse  data  from  an  embedded  sensor  with  probabilistic  life  predictions 
obtained  from  DARWIN.  These  fused  predictions  were  compared  with  measurements 
from  laboratory  demonstration  of  the  thin- film  sensor  under  the  DARPA  Material 
Prognosis  Program.  The  demonstration  consisted  of  a  growing  fatigue  crack  in  a  Ti-6A1- 
4Ti  beam  specimen  that  was  being  monitored  with  a  magnetostrictive,  thin- film 
(ultrasonic)  wireless  sensor  as  described  in  Ref.  [54], 

6.4.1  Crack  Detection  Experiment  with  the  Thin-Film  Sensor 

The  experimental  set-up  is  shown  in  Figure  106.  The  Ti-6A1-4V  three-point  bend 
specimen,  containing  a  central  surface  crack,  was  subjected  to  fatigue  loading.  Note  that 
the  sensor  was  placed  on  the  opposite  side  of  the  crack  to  demonstrate  that  bulk  waves 
were  being  generated.  Crack  size  measurements  were  made  at  both  25  °F  and  500  °F. 
Results  from  the  sensor  were  transmitted  both  wired  and  wirelessly.  The  sensor  was  also 
powered  through  a  wired  power  supply  since  energy  harvesting  has  not  yet  been  developed 
to  power  the  sensors.  The  surface  crack  was  generated  by  first  introducing  a  small  electro¬ 
discharge  machine  notch  followed  by  fatigue  pre-cracking  to  extend  the  crack.  During 
crack  extension,  marker  bands  were  periodically  generated  on  the  fracture  surface,  as 
shown  in  Figure  107,  by  increasing  the  mean  load  and  decreasing  the  load  range  to  alter 
the  crack  surface  morphology.  These  marker  bands  provided  precise  information  on  crack 
depth,  shape,  and  area  as  summarized  in  Table  10. 
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Figure  106.  Experimental  Setup  of  Fatigue-Cracked  Beam  Being  Monitored 

with  Thin-Film  Sensor 


Figure  107.  Fracture  Surface  of  Fatigue  Crack  in  Ti-6A1-4V 
Showing  Fatigue  Marker  Bands 
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Table  10.  Crack  Sizes  and  Areas  Measured  from  the  Surface  Crack  in  Figure  107 
That  Were  Detected  Using  the  Thin-Film  Sensor  Operating  at  500  °F 


Crack  Length  (in.) 

Crack  Depth  (in.) 

Crack  Area  (in.2) 

0.103 

0.033 

0.0034 

0.173 

0.036 

0.0062 

0.234 

0.039 

0.0091 

0.300 

0.05 

0.015 

0.420 

0.062 

0.026 

6.4.2  Thin-Film  Sensor  Uncertainty  Model 

The  signal-to-noise  ratio  obtained  from  the  wirelessly-transmitted,  TFS  results  at  25  °F  and 
500  °F  were  compared  against  known  crack  area  obtained  from  observation  of  post¬ 
inspection  marker  bands  on  the  crack  surface  (see  Table  10  and  Figure  107).  It  was 
determined  that  results  at  both  temperatures  were  similar,  thus  all  data  were  combined  to 
characterize  the  sensor  uncertainty  on  crack  size  determination.  Regression  results,  shown 
in  Figure  108,  were  used  to  quantify  the  uncertainty  in  the  sensor-derived  crack  area.  The 
regression  results  indicate  a  constant  standard  deviation  of  approximately  2349  mil2.  For 
analysis  purposes,  the  TFS  result  was  truncated  at  a  threshold  crack  size  of  2670  mil2 
(S/N  =  1),  i.e.,  no  size  information  is  provided  below  this  threshold  crack  size  since  for  this 
size  crack  the  sensor  signal  could  not  be  readily  distinguished  from  the  background  noise. 
This  threshold  size  value  was  set  at  -2  standard  deviations  which  set  the  mean  size  as 
7370  mil2  (for  a  =  c  =  66  mil).  In  summary,  the  crack  size  density  estimate  from  the  sensor 
response  was  modeled  as  a  normal  distribution  with  mean  7370  mil2  and  truncated  in  the 
left  tail  at  2670  mil2. 


Figure  108.  Thin-Film  Sensor  Experimental  Regression  Results  Used  to 
Establish  Sensor  Detection  Threshold  and  Measurement  Uncertainty 
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It  is  significant  to  note  that  the  minimum  detected  crack  area  of  2670  mils  ,  which  would, 
for  example,  correspond  to  a  semicircular  surface  crack  size  of  about  5 1  mils  (2  mm)  is 
2.5  times  smaller  in  size  and  15  times  smaller  in  area  than  the  crack  size  shown  to  be 
effective  in  maintaining  low  failure  rates  based  on  the  probabilistic  simulations  in 
Section  6.3  (Figure  103).  These  results  are  encouraging  since  they  suggest  that  the  current 
thin-film  embedded  sensing  technology  may  have  adequate  margin  to  enable  the  successful 
transition  from  the  lab  to  the  turbine  engine,  or  other  applications,  with  temperature  of 
500  °F  or  less. 

6.4.3  Initial  Crack  Size  Distribution 

A  central  issue  in  any  probabilistic  fracture  mechanics  analysis  is  the  initial  crack  size 
distribution.  Under  an  FAA-sponsored  effort,  the  turbine  engine  manufacturers  have 
developed  a  crack  size  distribution  for  machined  holes  in  disks  [42],  These  results  are 
shown  in  Figure  109,  which  represents  results  that  can  be  expected  for  a  fleet  of  disks.  The 
probability  density  is  “exponential-like”  in  that  it  has  a  high  probability  of  small  cracks 
with  a  long  right  tail  to  capture  the  remote  possibility  of  large  cracks.  This  distribution  is 
used  for  the  demonstration  problem.  As  will  be  shown,  the  initial  crack  size  distribution  is 
modified  significantly  based  upon  the  sensor  crack  size  estimation. 
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Figure  109.  Initial  Crack  Size  Distribution  for  Flaws  at  Machined  Holes  from  Ref.  [42] 
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6.4.4  Uncertainty  in  Fatigue  Crack  Growth  Rates 

The  uncertainty  in  fatigue  crack  growth  rates  was  captured  by  analyzing  an  extensive  data 
base  on  Ti-4A1-4V  obtained  in  an  Advisory  Group  for  Aerospace  Research  and 
Development  (AGARD)  round  robin  test  program  [55].  These  results  are  summarized  in 
Figure  1 10,  along  with  results  probability  density  functions  (PDFs)  at  various  fixed  values 
of  AK.  As  can  be  seen  in  the  figure,  the  CO  Vs  of  these  distributions  varied  as  a  function  of 
AK.  Specifically,  they  increased  at  both  low  and  high  AK  values.  However,  in  the 
intermediate  region  of  interest  in  the  current  simulation  the  COV  was  relatively  constant 
and  equal  to  20%  on  average.  A  value  of  50%  was  also  used  in  the  computations  to 
simulate  the  larger  uncertainty  that  would  be  associated  with  different  heats  of  material,  as 
would  be  the  case  in  fusing  data  in  an  actual  component. 


AK,  MPa/m 

Intermediate  growth  rate  regime  is  of  most  interest  here. 


Figure  110.  Fatigue  Crack  Growth  Data  in  Ti-6A1-4V  from  AGARD  Round  Robin 
Study  [55],  Along  with  PDFs  at  Various  AK  Values 
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6.4.5  Fused  Sensor  and  Fracture  Mechanics  Prediction  Results 


The  method  used  to  fuse  the  sensor  and  life  prediction  information  is  illustrated  in 
Figure  111.  An  sample  is  randomly  selected  from  the  initial  flaw  size  distribution 
(Figure  109)  and  appropriate  loading  cycles  (23,500  cycles)  are  applied  to  the  DARWIN 
life  prediction  software  to  grow  the  crack  to  the  point  of  the  initial  sensor  measurement  in 
the  experiment.  Since  the  sensor  predictions  at  this  point  are  in  terms  of  crack  area,  the 
crack  aspect-ratio  computed  by  DARWIN  is  fed  to  the  sensor  model  to  enable  conversion 
of  the  sensor-predicted  area  to  a  crack  depth  and  length.  This  procedure  is  performed 
numerous  times  using  Monte  Carlo  simulation  with  all  of  the  random  variables  indicated 
by  the  inner  loop  in  Figure  111.  The  simulation  results  are  then  used  to  formulate  a  non- 
parametric  PDF  of  crack  length  from  the  sensor  and  predictions  are  then  formulated  using 
the  adaptive  kernel  method  previously  described  in  Section  4.2.  The  above  procedure  is 
then  repeated  for  each  of  the  crack  lengths  where  sensor  measurements  were  made  in  the 
experiment  using  the  Bayesian  updating  procedure;  this  process  is  indicated  by  the  outer 
loop  in  Figure  111.  The  above  simulation  procedure  result  is  predicted  PDFs  as  shown  in 
the  lower  left  of  Figure  111  for  the  sensor  (blue),  life  prediction  model  (red),  as  well  as  for 
the  fused  results  (green). 
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Figure  111.  Method  Used  to  Fuse  Sensor  Measurements  and  DARWIN  Life  Predictions, 
and  with  Example  of  Resulting  Probability  Density  Functions 
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The  Bayesian  analysis  at  each  cycle  count  was  based  on  10,000  Monte  Carlo  simulations 
with  a  Standard  Deviation  of  0.0023  in.2  for  the  sensor  model,  and  COV  values  of  either 
20%  or  50%  for  life  scatter  -  that  was  assumed  to  be  controlled  by  the  fatigue  crack 
growth  rate  uncertainty.  The  results  for  the  simulations  with  the  life  scatter  of  20%  in  the 
DARWIN  predictions  are  shown  in  Figure  1 12  for  applied  cycles  of  23,500,  30,500,  and 
39,500  corresponding  to  the  cycle  counts  for  which  sensor  measurements  and  marker 
bands  were  applied  in  the  experiment.  As  shown  in  the  figure,  the  initial  Bayesian  results 
at  23,500  cycles  favor  the  DARWIN  predictions  because  of  the  relatively  low  COV  of  the 
DARWIN  predictions  relative  to  the  sensor  model  predictions.  However,  the  mean  values 
of  both  of  these  distributions  exceed  the  actual  surface  crack  length,  as  well  as  that 
predicted  by  the  sensor  model.  At  both  30,500  cycles  and  39,500  cycles,  the  DARWIN 
predictions  continue  to  over-predict  the  actual  crack  length  which  is  now  close  to  the 
predictions  of  the  sensor  model,  as  well  as  the  Bayesian  analysis.  This  over  prediction  by 
the  DARWIN  analysis  is  likely  due  to  the  fact  that  the  mean  of  the  fatigue  crack  growth 
rate  distribution  is  larger  than  that  of  the  actual  heat  of  Ti-6A1-4V  used  in  the  experiments. 
Thus,  the  availability  of  the  sensor  date  serves  to  in  this  example  resist  the  tendency  of  the 
DARWIN  analysis  to  over-predict  the  crack  length  and  consequently  the  Bayesian 
updating  provides  a  significantly  better  prediction  of  the  actual  crack  length  in  the 
experiment. 


Figure  112.  Comparison  of  Actual  Crack  Lengths  (Dashed  Lines)  with  Bayesian 
Predictions  (Green)  from  Fusion  of  Sensor  Model  Results  (Blue)  and  Fracture 
Mechanics  Model  Predictions  (Red)  with  Crack  Growth  COV  =  20% 
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The  results  for  the  simulations  where  the  life  scatter  was  increased  by  2.5  X  to  50%  (versus 
20%  in  the  analyses  shown  in  Figure  1 12)  are  shown  in  Figure  113.  The  overall  trends  in 
these  results  are  similar  to  those  for  the  smaller  life  scatter  case  in  Figure  1 12.  However, 
the  increased  uncertainty  in  the  DARWIN  analysis  eventually  results  in  all  of  the 
predictions  having  greater  uncertainty  as  reflected  in  the  broader  distributions. 
Nevertheless,  the  mean  values  of  the  sensor  and  Bayesian  predictions  remain  closer  to  the 
actual  crack  length  than  the  mean  value  of  the  DARWIN  analysis  as  the  result  of  the 
smaller  influence  of  the  DARWIN  analysis  because  of  its  higher  uncertainty. 

The  results  in  summarized  in  Figures  112  and  113  demonstrate  the  advantage  of  making 
predictions  having  both  sensor  information  and  analytical  life  predictions,  as  well  as  the 
benefit  of  using  Bayesian  analysis  in  fusing  both  sources  of  information. 


Figure  113.  Comparison  of  Actual  Crack  Lengths  (Dashed  Lines)  with  Bayesian 
Predictions  (Green)  from  Fusion  of  Sensor  Model  Results  (Blue)  and  Fracture 
Mechanics  Model  Predictions  (Red)  with  Crack  Growth  COY  =  50% 
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7. 


Conclusions 


The  following  points  summarize  the  key  accomplishments  and  significant  conclusions  that 

can  be  drawn  from  the  current  program: 

1 .  A  new  FaNG  model  was  developed  based  on  the  use  of  small  fatigue  crack  growth 
behavior  to  link  crack  nucleation  and  classical  large-crack  growth  regimes;  the  model 
is  particularly  useful  for  fatigue  analyses  involving  steep  stress  gradients  -  for 
example,  residual  stresses,  sharp  notches,  foreign  object  damage,  fretting  fatigue,  and 
thermal  stresses. 

2.  The  FaNG  Model  was  validated  by  comparing  predicted  versus  measured  life  of 
notched  fatigue  specimens,  and  peened  fatigue  specimens  of  Ti-6A1-4V  for  which 
residual  stress  measurements  were  also  available,  thereby  enabling  crack  growth  to 
be  computed  under  the  influence  of  the  steep  near-surface  residual  stress  field. 

3.  The  FaNG  Model  was  demonstrated  to  give  results  consistent  with  classical  analyses 
for  mild  stress  gradients. 

4.  Computing  fatigue  lives  in  a  fracture  critical  component  feature  (i.e.  bolt  hole  in  a 
generic  compressor  disc)  based  on  TACs  was  shown  to  be  conservative  by  a  factor  of 
two  to  five  compared  to  computed  fatigue  lives  based  on  actual  usage  from  a  variety 
of  mission  types.  These  results  provide  further  documentation  of  the  potential 
benefits  of  using  actual  mission  usage  instead  of  TACs  to  track  fatigue  damage, 
schedule  inspections  and  maintenance,  and  manage  the  life  fracture  critical 
components. 

5.  A  multi-component  da/dN(AK)  model  was  developed  and  used  to  quantify  fatigue 
crack  growth  rate  uncertainty  in  Ti-6A1-4V.  This  model  will  facilitate  the  explicit 
incorporation  of  fatigue  crack  growth  rate  material  properties  in  probabilistic 
analysis,  which  can  aid  in  materials  development,  as  well  as  materials  selection  for 
turbine  engine  applications. 

6.  Enhanced  probabilistic  methods  were  developed  that  can  increase  the  efficiency  and 
effectiveness  of  probabilistic  analyses  including: 

a.  Quantification  of  the  computational  benefits  of  parallel  processing  for 
probabilistic  fracture  mechanics  analysis  of  military  engines,  and 

b.  Development  and  validation  of  efficient  analytical  methods  for  computing 
probabilistic  sensitivity  factors  using  readily  available  computational  results. 

7.  The  utility  of  a  non-parametric  distribution  function,  based  on  an  adaptive  kernel 
method,  for  characterizing  uncertainty  in  mission  usage  data  was  demonstrated.  This 
method  can  be  used  to  automatically  update  mission  usage  as  additional  data  are 
captured  over  time. 
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8.  An  automated  method  for  PMI  was  developed  using  the  adaptive  kernel  method 
Initial  validation  of  the  PMI  method  for  automated  identification  of  mission  type 
from  given  usage  data  was  successfully  demonstrated  for  several  different  mission 
types.  This  tool  will  facilitate  the  automated  characterization  of  uncertainty  in 
mission  usage  and  thereby  lead  to  better  forecasting  of  the  impact  of  changes  in 
mission  mix  on  evolving  fatigue  damage. 

9.  Using  the  above  tool,  a  conceptual  design  was  formulated  to  synthesize  mission 
usage  and  forecast  the  effect  of  changes  in  mission  or  mission  mix  on  component 
reliability. 

10.  The  benefits  of  Bayesian  updating  for  fusing  information  from  on-board,  embedded 
sensing  and  probabilistic  life  predictions  were  demonstrated.  This  demonstration 
provided  further  evidence  for  the  feasibility  of  real-time,  or  near  real-time,  prognosis 
of  the  reliability  of  fracture  critical  turbine  engine  components. 

1 1 .  The  above  results  also  demonstrate  that  monitoring  component  fatigue  damage  (i.e. 
cracking)  with  embedded  sensors  having  measurement  sensitivities  that  are  five  times 
less  sensitive  than  those  used  for  depot  inspection  are  beneficial  for  ensuring  high 
reliabilities  in  fatigue  and  fracture  critical  turbine  engine  components  since  the  lower 
sensitivity  is  compensated  for  by  the  ability  of  embedded  sensors  to  continually 
interrogate  the  component  damage  state  -  for  example,  differential  damage 
measurements  acquired  for  each  flight. 

12.  Many  of  the  enhanced  life  prediction  and  probabilistic  analysis  methods  developed  in 
the  current  program  have  been  transitioned  to  the  commercial  computer  code 
DARWIN,  which  is  available  royalty-free  to  the  government  and  through  commercial 
licensing  to  industry;  these  methods  which  are  accessible  through  the  user  friendly 
GUI  include: 

a.  FaNG  model 

b.  Residual  stress  analysis  for  surface  cracks 

c.  Enhanced  mission  mixing 

d.  Probabilistic  sensitivities. 

13.  The  current  work  identified  several  technical  challenges  to  hybridizing  physics-based 
probabilistic  analysis  (DARWIN)  and  data-driven  (ProDAPS)  including: 

a.  Physics-based  models  are  available  for  fracture-critical  components,  while 
data  driven  methods  have  not  been  developed  and  matured  for  the  same 
damage  modes 

b.  For  those  limited  cases  where  the  two  methods  are  compatible,  data  are  often 
lacking  because  of  inaccessibility  to  software  needed  to  interpret  raw  data, 
incompatibility  between  different  data  bases,  and/or  a  general  lack  of 
infrastructure  for  storing  and  processing  large  amounts  of  data. 
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8. 


Recommendations 


1 .  Explore  the  use  of  the  probabilistic  version  of  the  FaNG  model  to  elucidate  the 
bimodal  distribution  of  coupon  S-N  fatigue  data.  This  will  require  a  consistent  set  of 
representative  material  properties  [S-N,  da/dN(AK),  and  small  crack  data  (optional)], 
as  well  as  initial  defect  distributions  inherent  to  the  material. 

2.  Perform  further  validation  of  the  FaNG  model  over  a  broader  range  of  materials  (e.g. 
Ni-base  alloys),  stress  gradients  (geometries,  thermal  stresses,  and/or  beneficial 
residual  stresses).  This  will  require  representative  material  properties 

[S-N],  as  well  as  coupon,  subcomponent  or  component  test  data  with  known  local 
stress  distributions  at  fatigue  critical  locations,  known  loading  history,  as  well  as 
characterized  residual  stresses  -  if  they  exist  in  the  coupon  or  component. 

3.  Perform  further  assessment  of  the  viability  of  the  PMI  to  accurately  identify  mission 
type  for  a  given  usage  (time  history);  testing  robustness  of  the  method  will  require 
access  to  more  usage  data  from  a  wider  range  of  mission  types  for  which  the  mission 
type  is  known. 

4.  Add  the  multi-component  fatigue  crack  growth  model  to  DARWIN  so  that  the 
software  code  can  be  used  as  a  research  tool  to  examine  the  potential  impact  of  new 
materials,  with  improved  crack  growth  behavior,  on  turbine  engine  reliability. 

5.  Implement  current  Air  Force  plans  to  develop  data  storage  and  data  mining  tools,  as 
well  as  the  ability  of  share  data  among  the  existing  Air  Force  data  bases  for  engines, 
as  well  as  airframes.  This  would  facilitate  the  validation  of  existing  and  new 
prognostic  and  diagnostic  tools.  It  is  also  important  to  develop  procedures  to  allow 
researchers  to  access  data  while  preserving  the  confidentiality  of  the  information. 
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OPS 
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Probabilistic  Density  Function 
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Probability  of  Detection 
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Power  Lever  Angles 
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Probabilistic  Mission  Identification 
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Program  Research  and  Development  Announcement 
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Probabilistic  Diagnostic  and  Prognostic  System 
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SC 
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SP 

SwRI 
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TACs 
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UTSA 
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Q 

E 

Fj 

K 

Kc 

KIC 


Revolutions  Per  Minute 
Residual  Stress 

Society  of  Automotive  Engineering 
Small  Crack 

stress  versus  number  of  cycles  to  failure 
Shot-Peening  (a  surface  treatment) 

Southwest  Research  Institute 
Smith- W  atson-Topper 
Total  Accumulated  Cycle(s) 

Thin-Film  Sensor 
Turbine  Rotor  Material  Design 
United  States  Air  Force 
University  of  Texas  at  San  Antonio 

Description 

Pre-exponential  constants  in  multicomponent  crack  growth  law 
Young’s  modulus,  ksi 
jth  fleet 

Applied  stress  intensity  factor  range 

Fracture  toughness  measured  during  a  fatigue  crack  growth  test  under 
non-plane  strain  conditions  (typically  Kc  >  Kjc) 

Plane-strain  fracture  toughness 
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Mi 

Mi 

Nlc 

Nsc 

Nf 

Nf* 

Ni 

NiSWT 

Ntn 

Pf 

R 

S/N 

a 

ao 

da/dN 

/ M+G 

ri+1 

Jm 

fs 

h 


Exponents  in  multicomponent  crack  growth  law 
Cycles  of  large  crack  growth 
Cycles  of  small  crack  growth 
Cycles  to  failure 

Flights  to  failure  at  100%  Design  RPM 
TACs  per  Flight 

Cycles  to  initiation 

Cycles  to  initiation  computed  from  Smith- Watson-Topper  criterion 
Cycles  to  “true”  crack  nucleation 
Probability  of  Failure 

Load  (or  stress)  ratio  used  to  characterize  fatigue  loading 
Signal-to-noise  ratio 
Confidence  bound 
Crack  size 

Small  fatigue  crack  parameter 
Fatigue  crack  growth  rate,  in./cycle 

Bayesian  updated  PDF  for  next  flight  (input  to  model) 

Model  PDF  for  i+1  growth  cycle 

Sensor  PDF  for  flight  i 

Kernel  window  width 
Number  of  data  points 
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A K 
A Keq 

A Kth 
A  ae 
O 

£(i,max 

(7 

Cm  ax 

z 

80 

O0Z 

80 


Applied  stress  intensity  factor  range 

Equivalent  stress  intensity  factor  range  that  accounts  for  increased  crack¬ 
driving  force  of  small  fatigue  cracks 

Fatigue  threshold  stress  intensity  factor  range,  ksiVin. 

Fatigue  endurance  limit  in  a  smooth  specimen 

Standard  normal  cumulative  distribution  function 

Maximum  principal  strain  amplitude 

Mean  value  of  a  distribution 

Standard  deviation  of  a  distribution 

Maximum  normal  stress  on  the  critical  plane  of  ea,max 

Probabilistic  sensitivity  of  random  variable  0  with  respect  to  uncertainty 
in  mean  crof  the  response  Z 

Probabilistic  sensitivity  of  random  variable  0  with  respect  to  uncertainty 
in  mean  //  of  the  response  Z 
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